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Abstract 

The  US  Air  Force’s  ability  to  protect  space  assets  is  enhanced  by  a  proficiency 
in  satellite  proximity  operations  and  Space  Situational  Awareness  (SSA).  In  pursuit 
of  that  proficiency,  this  research  develops  a  key  capability  of  interest  to  mission  plan¬ 
ners;  the  ability  of  a  deputy  satellite  to  “hover”  within  a  defined  volume  fixed  in 
the  vicinity  of  a  chief  satellite  for  an  extended  period  of  time.  Previous  research  has 
developed  initial  methodologies  for  maintaining  restricted  teardrop  hover  orbits  that 
exist  in  a  plane  fixed  within  the  chief’s  local  reference  frame.  These  methods  use 
the  natural  drift  of  the  deputy  satellite  in  the  relative  frame  and  impulsive  thrust 
to  keep  the  deputy  in  a  bounded  volume  relative  to  the  chief,  but  do  not  address 
fuel-optimality.  This  research  extends  and  enhances  that  work  by  finding  optimal 
trajectories,  produced  with  discrete-thrusts,  that  minimize  fuel  spent  per  unit  time 
and  stay  within  the  user-defined  volume,  thus  providing  a  practical  hover  capability  in 
the  vicinity  of  the  chief.  The  work  assumes  the  Clohessy-Wiltshire  closeness  assump¬ 
tion  between  the  deputy  and  chief  is  valid,  however,  elliptical  chief  orbits  are  allowed. 
Using  the  new  methodology  developed  in  this  work,  feasible  closed  and  non-closed 
relative  orbits  are  found  and  evaluated  based  on  a  fuel  criterion  and  compared  to 
an  easily  calculated  continuous-thrust  baseline.  It  is  shown  that  in  certain  scenarios 
the  discrete-thrust  solution  provides  the  lowest  overall  fuel  cost.  These  scenarios  are 
generally  constrained  to  a  smaller  total  time-of- flight.  A  simple  check  is  proposed 
that  enables  the  mission  planner  to  make  the  correct  strategy  choice. 
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Optimal  Control  Strategies  for 
Constrained  Relative  Orbits 

I.  Introduction 

The  dawn  of  the  20th  century  saw  humanity’s  hrst  hesitant  steps  from  the 
surface  of  Earth,  setting  off  a  technological  and  engineering  explosion  that  put  a 
manmade  object  into  space  a  mere  54  years  after  that  first  flight.  The  US  Military 
has  made  enormous  strides  in  utilizing  the  space  environment  to  provide  capability 
and  multiply  overall  military  effectiveness  in  combat  operations.  As  the  world’s  tech¬ 
nological  base  has  grown,  our  allies  and  enemies  are  also  taking  advantage  of  the 
opportunities  gained  with  their  own  assets  in  space.  Although  doctrine  is  changing, 
space  has  historically  been  viewed  mostly  as  a  medium  of  non-interference  due  in  large 
part  to  an  inability  to  affect  assets  on  orbit.  This  could  be  viewed  as  a  modification 
of  the  “Big  Sky”  theory  of  early  powered  flight  in  which  opposing  forces  merely  waved 
at  each  other  as  they  passed  by  on  their  way  to  the  fight.  In  today’s  military  envi¬ 
ronment  this  simple  philosophy  is  no  longer  viable.  Recent  events  have  proven  that 
our  space  assets  are  not  beyond  our  enemy’s  reach  and  therefore  must  be  protected. 
This  new  mission  for  the  USAF  is  called  Counterspace  and  is  defined  in  Air  Force 
Doctrine  Document  1  [10]  as 

those  kinetic  and  nonkinetic  operations  conducted  to  attain  and  maintain 
a  desired  degree  of  space  superiority  by  the  destruction,  degradation,  or 
disruptions  of  enemy  space  capability 

Counterspace  is  separated  into  two  pieces:  Offensive  Counterspace  (OCS)  and 
Defensive  Counterspace  (DCS).  DCS  preserves  space  capabilities  from  enemy  threats 
while  OCS  operations  seek  to  affect  non-US  space  assets  negatively.  On  the  DCS 
side,  there  are  a  number  of  scenarios  in  which  it  would  be  advantageous  for  a  friendly 
micro-satellite  to  stay  within  the  local  area  of  a  larger  friendly  satellite  in  a  protection 
or  inspection  role.  The  feasibility  of  this  type  of  mission  using  closed  orbits  with  the 
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target  satellite  in  the  center  was  investigated  by  Thomas  et  al  [54],  The  protective 
“sentinel”  mission  may  require  the  microsat  to  remain  in  an  orbit  near  the  target 
satellite  but  transition  to  a  dehned  relative  location  for  defensive  mode  upon  indication 
of  attack. 

In  addition  to  DCS  measures,  there  are  a  host  of  capabilities  that  will  be  essential 
for  maintaining  critical  space  systems  in  the  future.  These  include  repairing,  refueling, 
upgrading,  augmenting,  and  otherwise  servicing  on-orbit  space  assets.  The  linchpin 
for  all  of  these  capabilities  is  the  ability  to  perform  close-proximity  operations.  A 
specific  type  of  close-proximity  operation,  as  mentioned  briefly  above,  is  the  ability 
to  “hover”  in  a  specific  zone  relative  to  the  target  satellite.  Hover  capability  has  been 
demonstrated  for  a  few  constrained  cases  (i.e.,  in  specific  places  relative  to  the  target 
satellite)  or  missions  about  smaller  planetary  bodies  such  as  asteroids.  The  goal  of 
this  dissertation  research  is  to  extend  hover  capability  to  anywhere  within  the  target 
satellite’s  local  area  and  for  a  target  satellite  on  any  closed  Keplerian  orbit.  Note  that, 
since  hovering  will  typically  not  be  on  a  natural,  drift-free  relative  orbit,  thrusting 
will  be  required.  Thus  the  research  will  focus  on  maximizing  the  time  of  hover  for 
a  given  amount  of  fuel  and,  in  doing  so,  will  address  the  feasibility  of  a  variety  of 
potential  scenarios  under  consideration  by  military  planners.  For  the  purposes  of  this 
research,  the  target  satellite  of  interest,  which  is  located  at  the  center  of  the  relative 
frame,  will  be  called  the  “chief’  satellite  and  the  satellite  operating  in  proximity  to 
the  chief  is  the  “deputy” .  The  problem  statement  addressed  in  this  dissertation  is  to 

develop  a  control  strategy  to  place  a  deputy  satellite  inside  a  specific  lobe 
defined  in  the  chief  body-fixed  frame  and  keep  it  there  in  a  fuel-optimal 
manner. 

Solving  this  problem  will  allow  us  to  answer  two  questions: 

1.  Can  a  discrete-thrust  trajectory  be  found  that  outperforms  a  benchmark  continuous- 
thrust  solution  for  a  fuel  criterion  of  optimality? 

2.  Can  we  quickly  and  robustly  estimate,  with  reasonable  accuracy,  the  amount  of 
AV  required  to  stay  within  a  specific  lobe? 
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The  following  chapters  detail  a  solution  method  and  analyze  results  in  order 
to  answer  those  questions.  Chapter  II  presents  an  overview  of  proximity  operation 
and  relative  motion  research  currently  in  the  literature.  Chapter  III  introduces  the 
elements  needed  to  solve  the  problem  including  the  equations  of  motion,  development 
of  cost  and  constraint  functions,  and  analysis  of  symmetries  in  the  problem  that 
reduce  the  number  of  cases  needed  to  confirm  the  conclusions.  Chapter  IV  describes 
how  the  optimal  trajectory  is  found,  along  with  definitions  of  the  initial  and  final 
conditions.  This  chapter  ends  with  an  overview  of  the  research  cases  found  in  the 
results.  Chapter  V  presents  and  discusses  optimal  trajectories  for  four  different  lobes. 
Finally,  conclusions  and  suggestions  for  future  research  are  found  in  Chapter  VI. 
The  appendices  are  ordered  such  that  they  build  from  mathematical  preliminaries  to 
derivation  of  the  foundation  equations  to  application  of  those  equations.  Therefore, 
they  will  not  be  referenced  in  alphabetical  order  in  the  main  document. 
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II.  Literature  Review 


2.1  Background 

It  was  G.W.  Hill  who  developed  the  restricted  three-body  problem  of  the  Sun- 
Earth-Moon  system  in  terms  of  a  relative  rotating  reference  frame  [16].  The  advent 
of  artificial  satellites  and  the  potential  for  constructing  larger  structures  in  space 
and/or  docking  with  those  structures  necessitated  a  method  of  rendezvous.  This 
need  led  to  W.H.  Clohessy  and  R.S.  Wiltshire’s  adoption  of  Hill’s  methodology  and 
the  production  of  their  famed  Clohessy-Wiltshire  relative  equations  of  motion  in  the 
early  1960’s  [8].  These  equations  allow  for  not  only  docking  and  rendezvous  but  also 
close-proximity  operations  in  which  a  deputy  satellite  is  placed  in  a  closed  relative 
orbit  about  the  chief.  Although  the  Clohessy-Wiltshire  equations  (abbreviated  CW 
henceforth)  are  valid  only  for  chief  satellites  in  a  circular  orbit,  they  have  proven 
quite  useful  in  solving  a  wide  variety  of  satellite  dynamics  problems.  In  addition, 
by  assuming  the  chief  and  deputy  are  close  in  comparison  to  the  orbital  radius  of 
the  chief,  the  CW  equations  very  nicely  linearize  to  a  form  that  is  suitable  for  linear 
analysis  and  control  techniques.  These  restrictions  can  be  removed  with  more  complex 
sets  of  equations  and  have  been  studied  extensively  in  recent  years  [5,20,25]. 

Research  into  relative  orbit  dynamics  (see  references  in  Section  2.2  to  Sec¬ 
tion  2.6)  has  exploded  over  the  last  two  decades  as  the  potential  payoff  for  coor¬ 
dinated  satellite  formations  has  been  brought  into  sharp  focus.  Covernment  interest, 
specifically  from  the  USAF  and  DARPA,  is  clear  by  the  numerous  programs  dedi¬ 
cated  to  formation  flying  and  its  associated  technology  development.  The  TechSat 
21  program  [31],  although  now  defunct,  investigated  technologies  critical  to  satel¬ 
lite  formations  such  as  micro-satellite  bus  and  micro-propulsion.  TechSat  21  had 
planned  to  be  a  proof  of  concept  for  distributed  mission  architecture,  sparse  aperture 
sensing,  and  collaborative  behavior  [6,49].  In  an  effort  to  demonstrate  proximity  op¬ 
erations,  the  AFRL  Space  Vehicles  directorate  has  executed  the  XSS-10  and  XSS-11 
missions  [52].  The  XSS-11,  launched  in  April  2005,  successfully  performed  a  variety 
of  rendezvous  and  proximity  missions  of  several  US-owned,  dead,  or  inactive  space 
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objects.  DARPA’s  Orbital  Express  Space  Operations  Architecture  similarly  sought 
to  validate  a  host  of  proximity  operations,  including  autonomous  on-orbit  refueling 
and  reconfiguration  of  satellites  [40].  It  successfully  launched  on  8  March  2007  and 
completed  its  technology  demonstration  mission  on  29  June  2007.  Finally,  NASA’s 
Demonstration  of  Autonomous  Rendezvous  Technology  (DART)  program  launched 
in  April  2005  and  successfully  completed  the  location  and  rendezvous  phases  of  its 
operations  but  was  unable  to  complete  all  of  its  close-proximity  and  circumnavigation 
operations  due  to  lack  of  fuel  [1,39]. 

The  technical  papers  have  been  collected  in  five  categories: 

1.  Dynamic  Analysis 

2.  Formation  Establishment,  Maintenance,  and  Reconfiguration 

3.  Effect  of  Perturbations 

4.  Extension  to  Elliptical  Chief  Orbits 

5.  Constrained  Formation  Geometries 

2.2  Dynamic  Analysis 

The  first  group  consists  of  papers  that  investigated  either  alternate  means  of 
expressing  the  relative  equations  of  motion  or  higher  order  expansions  of  the  CW 
equations.  Many  researchers  have  abandoned  the  relative  Cartesian  coordinates  of 
the  CW  equations  for  those  based  on  classical  orbital  elements  of  the  chief  or  differ¬ 
ences  between  the  chief  and  deputy  elements  [45,46].  This  has  been  used  to  derive 
minimum,  maximum,  and  mean  distance  between  the  two  spacecraft  [15].  An  approx¬ 
imate  second-order  solution  to  the  relative  orbit  equations  was  derived  by  Karlgaard 
and  Lutze  [22]  that  shows  a  two-orders-of-magnitude  improvement  over  the  linear 
solution  over  one  period  of  the  reference  orbit.  A  third-order  analytical  solution  was 
developed  by  Richardson  and  Mitchell  [43].  In  an  effort  to  compare  the  various  rel¬ 
ative  orbit  models  head  to  head,  Alfriend  and  Yan  [2]  have  created  a  useful  tool  for 


5 


evaluating  and  comparing  the  accuracy  of  different  models  through  the  use  of  a  mod¬ 
eling  error  index.  Finally,  Amico  and  Montenbruck  [9]  have  adapted  the  concept  of 
eccentricity /inclination- vector  separation  (originally  used  for  geostationary  satellites) 
to  low  Earth  orbit  formations  and  show  its  utility  for  proximity  analysis  as  well  as 
orbit  control. 

2.3  Formation  Establishment,  Maintenance,  and  Reconfiguration 

Armed  with  the  right  set  of  equations,  researchers  have  looked  at  numerous 
techniques  for  placing  satellites  in  specific  formations  and  keeping  them  there.  Of 
course,  one  of  the  main  advantages  of  a  satellite  formation  is  its  ability  to  adapt  to 
new  missions,  upgrade  with  new  satellites,  and  gracefully  handle  the  loss  of  a  dam¬ 
aged  or  dead  member,  resulting  in  investigation  of  methods  of  time-  and  fuel-optimal 
formation  reconfigurations.  Vaddi  et  al.  [58]  derived  an  analytical,  two-impulse  solu¬ 
tion  using  orbital  element  differences  to  establish  and  reconfigure  a  circular  satellite 
formation.  Yeh  et  al.  [68]  used  sliding  mode  control  with  the  nonlinear  Hill’s  equa¬ 
tions  to  maintain  the  formation  in  the  presence  of  drag,  third-body  effects,  and  an 
oblate  Earth.  A  similar  approach  was  used  by  Massey  and  Shtessel  [32].  Lovell  and 
Tragesser  [30]  proposed  a  multiple-impulsive  maneuver  for  reconfiguration  based  on 
an  alternative  parameterization  of  the  CW  equations  which  allowed  mission  planners 
to  search  for  optimal  solutions.  Paiewonsky  and  Woodrow  [41]  tackled  the  problem 
of  finding  time-optimal  solutions  to  rendezvous  with  constraints  on  the  thrust  magni¬ 
tude  and  fuel  availability.  Guibout  and  Scheeres  [14]  used  a  Hamiltonian  approach  to 
solve  the  two-point  boundary  value  problem  of  formation  reconfiguration.  Palmer  [42] 
found  analytic  solutions  for  minimum-fuel  transfer  paths  between  two  relative  orbits 
with  a  fixed  time  of  flight  and  boundary  conditions.  Although  much  of  the  recent 
literature  has  focused  on  nonlinear  control  techniques,  the  original  work  and  some 
current  authors  utilize  linear  control  methods  [21,60].  Other  similar  works  on  this 
topic  are  found  in  [11,34]. 
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2-4  Effect  of  Perturbations 

The  majority  of  the  research  community  has  focused  on  oblate  Earth  and  at¬ 
mospheric  drag  perturbations  as  the  primary  influences  on  relative  satellite  orbits. 
Humi  and  Carter  [19]  investigated  relative  motion  in  the  presence  of  linear  drag.  This 
work  is  later  expanded  to  allow  for  quadratic  drag  [7].  Schweighart  and  Sedwick  [47] 
developed  a  set  of  linearized  constant  coefficient  differential  equations  that  describe 
satellite  relative  motion  in  the  presence  of  the  J2  potential.  They  later  updated  and 
simplified  the  cross-track  equations  in  reference  [48] .  Wiesel  [64]  employs  Floquet  the¬ 
ory  to  include  all  zonal  harmonics  of  the  Earth’s  gravitational  field,  which  produces 
accuracies  over  two  orders  of  magnitude  better  than  CW,  and  then  finds  an  optimal 
impulsive  control  law  for  stationkeeping  [65].  Leonard  et  al.  [27]  found  methods  of  us¬ 
ing  differential  drag  to  control  relative  positions.  Vadali  et  al.  [55]  determined  initial 
conditions  to  match  mean  J2-induced  drift  rates,  thereby  requiring  less  fuel  to  main¬ 
tain  the  formation.  A  second  method  imposed  periodic  boundary  conditions  on  the 
relative  position  and  velocity  in  a  rotating  coordinate  system.  Finally,  Williams  and 
Wang  [66]  looked  at  the  effects  of  solar  radiation  pressure  on  formations  in  highly  el¬ 
liptical  orbits.  Other  similar  works  on  dealing  with  perturbations  are  found  in  [35,56]. 

2.5  Extension  to  Elliptical  Chief  Orbits 

More  recent  research  has  investigated  control  of  the  nonlinear  equations  of  mo¬ 
tion  either  through  higher  order  approximations  or  attempts  to  handle  the  full  nonlin¬ 
ear  equations.  This  is  motivated  by  the  rather  restrictive  assumption  of  the  original 
CW  equations  that  the  chief  exists  on  a  circular  inertial  orbit.  Inalhan  et  al  [20] 
reintroduced  the  community  to  previous  extensions  of  the  CW  equations  to  eccentric 
chief  orbits.  They  also  developed  an  algorithm  to  find  initial  conditions  that  pro¬ 
duced  stable  periodic  solutions.  Yamanaka  and  Ankersen  [67]  found  a  set  of  linear 
differential  equations  with  time-dependent  coefficients  that  describe  relative  motion 
of  satellites  in  elliptical  orbits.  Alfriend  et  al.  [3]  found  a  second-order  theory  for 
relative  motion  that  allows  for  any  eccentricity  and  contains  first-order  J2  effects  that 
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can  be  easily  modified  for  higher  order  geopotential  terms.  Gim  and  Alfriend  [13] 
continue  this  exploration  by  finding  a  closed-form  state  transition  matrix  using  both 
osculating  and  mean  elements  that  allows  for  elliptical  orbits  and  J2.  Mitchell  and 
Richardson  [36]  have  developed  a  method  to  control  a  deputy  satellite  using  a  first- 
order  approximation  of  the  nonlinear  CW  equations  by  forcing  the  nonlinear  system 
onto  an  invariant  manifold  of  the  linear  system.  Along-track  position  drift  was  re¬ 
duced  by  a  factor  of  12  with  only  a  slight  control  increase  in  the  radial  and  cross-track 
drift.  Gurfil  and  Kasdin  [23]  present  a  method  to  obtain  higher-order  approximations 
of  the  relative  motion  in  which  the  coefficients  of  the  time  series  are  functions  of  the 
orbital  elements.  Broucke  [5]  is  the  first  to  present  results  with  time  as  the  indepen¬ 
dent  variable,  leading  to  a  set  of  linearized  equations  for  deputy  motion  close  to  a 
target  in  an  elliptic  orbit.  Lane  and  Axelrad  [25]  have  developed  a  set  of  geometric 
relationships  that  approximate  relative  motion  for  satellites  in  eccentric  orbits.  Vaddi 
et  al.  [57]  established  corrections  to  the  initial  conditions  of  the  linearized  equations 
of  motion  that  produce  periodic  solutions  for  the  nonlinear  GW  equations. 

2.6  Constrained  Formation  Geometries 

The  final  area  of  research  reviewed  is  an  emerging  class  of  problems  and  is  the 
basis  of  this  dissertation.  The  work  discussed  above  was  almost  exclusively  concerned 
with  relative  orbits  that  followed  closed  elliptical  paths  about  the  chief  satellite  and 
how  to  maintain  or  reconfigure  those  orbits.  As  discussed  in  the  introduction,  there 
may  be  cases  in  which  we  desire  the  deputy  not  to  orbit  around  the  chief  but  to  stay 
in  a  specific  constrained  volume  relative  to  the  chief.  Since  the  natural  dynamics  of 
the  system  indicate  zero-energy  closed-path  relative  orbits  exist  only  in  a  restricted 
case,  staying  within  an  arbitrary  constrained  volume  will  require  additional  energy 
(i.e.,  thrusting).  The  first  foray  into  finding  “hovering”  orbits  considered  finding 
closed  orbits  in  the  inertial  plane  of  the  chief  satellite.  Hope  and  Trask  [17]  found 
that  utilizing  the  natural  drift  of  the  relative  orbit  led  to  a  trajectory  in  the  inertial 
plane  that  intersected  itself,  providing  a  point  at  which  to  perform  a  single  impulsive 


thrust  to  keep  the  deputy  in  a  constrained  area.  This  work  was  expanded  by  Lovell 
and  Tollefson  [29]  who  developed  a  simpler  closed-form  solution  for  designing  the  size 
and  shape  of  the  hovering  orbit.  Finally,  there  has  been  work  in  optimal  control  of 
hovering  satellites  within  weak  gravitational  fields  such  as  missions  near  an  asteroid 
[4,12,18,44].  While  interesting,  the  low  gravity  field  assumption  has  little  applicability 
to  this  research. 

2. 7  Summary 

The  literature  is  replete  with  relative  satellite  motion  research  but  is  mostly 
concerned  with  higher  fidelity  modeling  of  the  equations  of  motion  or  creating  and 
maintaining  fixed  geometry  formations.  Initial  research  has  been  done  on  constrained 
relative  orbits  but  has  approached  the  problem  by  analyzing  trajectory  geometry  to 
find  feasible  orbits  without  addressing  the  fuel-optimality  of  those  orbits.  This  hereto¬ 
fore  unexplored  area  of  relative  satellite  motion  research  provides  ample  opportunity 
to  contribute  to  the  community. 
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III.  Methodology 


This  chapter  provides  the  mathematical  underpinnings  of  the  problem  and  pro¬ 
posed  solution.  Starting  with  the  equations  of  motion,  suitable  cost  and  constraint 
functions  are  derived.  Also  provided  is  a  lengthy  discussion  on  symmetries  within  the 
problem  that  greatly  reduce  the  number  of  examples  needed  to  confirm  the  results 
and  make  conclusions. 

Previous  work  on  “hovering”  orbits  has  considered  closed  orbits  in  the  inertial 
plane  of  the  chief  satellite  [17].  The  teardrop  orbit  is  designed  by  finding  a  drifting 
relative  orbit  in  the  inertial  plane  that  intersects  itself,  providing  a  point  at  which 
to  perform  a  single  impulsive  thrust  to  keep  the  deputy  in  a  constrained  area.  This 
research  adds  fuel-optimality  as  a  condition  while  finding  constrained  relative  orbits 
and  extends  the  problem  to  three  dimensions.  In  short,  the  problem  statement,  as 
introduced  in  Chapter  I,  is  to 

develop  a  control  strategy  to  place  a  deputy  satellite  inside  a  specific  lobe 
defined  in  the  chief  body-fixed  frame  and  keep  it  there  in  a  fuel-optimal 
manner. 

The  control  strategy  proposed  is  to  execute  impulsive  thrusts  such  that  the  location, 
magnitude,  and  direction  of  thrust  are  results  of  an  optimization  algorithm  developed 
herein.  The  simplest  realizable  case  is  the  one  in  which  two  satellites  share  the  same 
circular  orbit  but  have  different  angular  positions  in  their  inertial  orbits  about  the 
Earth.  In  the  absence  of  perturbations,  the  two  would  stay  in  a  fixed  relative  position 
with  each  other.  We  desire  more  flexibility  in  the  relative  placement,  however,  and 
seek  to  define  a  general  closed  lobe  (region)  of  arbitrary  size,  location,  and  orientation 
near  the  chief  satellite  and  fixed  in  the  relative  frame  that  bounds  the  relative  motion 
of  the  deputy  (see  Figure  1).  This  is  the  essence  of  “hovering”  and  is  formally  defined 
as 


Hovering:  Remaining  inside  a  specified  volume  defined  in  a  chief-centered 
reference  frame. 

We  especially  want  hovering  trajectories  that  are  fuel-optimal  which  we  define  as 
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Fuel-Optimal  Hovering  Trajectory:  A  trajectory  that  minimizes  the 
total  amount  of  fuel  spent  per  unit  time  of  hovering. 

We  also  want  to  evaluate  closed  and  non-closed  trajectories  where 

Closed  Trajectory:  A  trajectory  in  which  the  final  relative  position  and 
velocity  are  equal  to  the  initial,  thus  allowing  a  repeatable  relative  orbit. 

Figure  1  illustrates  a  general  lobe  with  a  center  that  is  defined  by  the  angles  a  and 
(3  and  a  distance  7  from  the  chief  satellite.  The  reference  frame  centered  on  the  chief 
satellite  is  defined  as  follows.  The  X  is  oriented  along  a  line  from  the  center  of  the 
Earth  to  the  chief,  Z  is  perpendicular  to  the  orbit  plane  of  the  chief  and  Y  completes 
the  frame  as  the  cross  product  Z  x  X.  The  in-track  direction  is  aligned  with  the 
velocity  vector  of  the  chief  when  in  a  circular  orbit.  This  frame  is  commonly  referred 
to  as  the  Local- Vertical/Local-Horizon  (LVLH)  frame. 


Figure  1:  General  Problem  Formulation 


3. 1  Notation 

The  following  notation  is  used  in  this  research.  The  positions,  Xj,  Ui,  and  Zi 
are  the  coordinates  of  the  ith  point  defined  in  the  LVLH  frame  and  indicate  where 
an  impulsive  thrust  may  occur.  An  ellipse  (see  Figure  2(a))  is  used  to  define  the 
lobe  boundary  in  the  chief’s  orbit  plane  {XY).  This  provides  significant  flexibility  for 
designing  hover  regions  without  overly  complicating  the  mathematics.  As  we  will  see 


11 


in  later  sections,  motion  in  the  Z  direction  is  decoupled  from  motion  in  the  XY  plane, 
thus  an  elliptical  cylinder  (see  Figure  2(b))  is  selected  in  order  to  prevent  recoupling 
of  those  equations.  The  lobe  center  is  located  at 


xt  =  ^xY  0L  =  ^  cos  a  sin  P 
Ul  =  JxY  =  7  sin  a  sin  P 

zl  =  7  cos  P 


where  a  is  the  angular  position  of  the  lobe  center  in  the  XY  plane  measured  counter¬ 
clockwise  from  the  X  axis,  P  is  the  angular  position  of  the  lobe  center  measured  from 
the  Z  axis,  7  is  the  distance  from  the  chief  to  the  lobe  center,  7^,^  is  the  projection 
of  7  in  the  XY  plane,  r]  is  the  angular  orientation  of  the  ellipse  measured  counter¬ 
clockwise  from  the  X  axis,  and  and  Ty  are  the  ellipse  axes.  For  lobes  that  exist  in 
three  dimensions  we  must  also  define  h,  the  half  height  of  the  elliptical  cylinder. 


Figure  2:  Lobe  Parameters 


Polar  coordinates  prove  useful  for  defining  impulsive  thrust  locations  on  the 
ellipse,  thus  the  angle  and  radius  of  the  zth  point  are  designated  by  pi  and  r*  where 
(derivation  in  Appendix  B) 

Vi  —  "f  sin  a  sin  P 
Xi  —  ^  cos  a  sin  P 

Ti  =  A/(a:j  —  ycosasin/?)^ -h  (r/j  —  ysinacos/d)^ 


Pi  =  tan  ^ 
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Conversion  back  to  cartesian  coordinates  yields 


Xi  =  ^  cos  a  sin  (3  + 
Vi  =  1  sin  a  sin  [3  + 


T^Ty  cos 


cos2  {'tpi  -ri)+T^  sin^  {'tpi  -  rj) 
_ TxTy  sin  xpi _ 

^ r2  cos2  -ri)  +  T^  sin^  -  rj) 


The  time-of-flight  between  the  ith  and  [i  +  l]th  points  is  denoted  by 


Note  that,  unless  otherwise  specihed,  the  deputy  starts  at  the  initial  position  at  t  =  0. 
Time-of-flight  is  easily  scaled  by  converting  it  to  fractions  of  the  chief’s  orbital  period 
(P)  defined  as  [62] 


P  =  2n.i^=—  (1) 

]j  H  n 

where  osma  is  the  semi-major  axis  of  the  chief  orbit,  fi  is  the  gravitational  parameter 
(constant  for  a  specific  two-body  system),  and  n  is  the  mean  motion  of  the  chief 


[I 

n 


G{mi  -I-  m2) 


^SMA 


where  G  is  the  fundamental  gravitational  constant  and  equal  to  6.6695  x  10 
and  m  is  mass.  Time-of-flight  as  a  fraction  of  chief  orbit  period  (T)  is 


P 


27r 


(2) 


Two  velocities  are  associated  with  each  impulsive  thrust  point,  an  arriving  velocity 
which  is  a  function  of  the  previous  position,  the  current  position,  and  the  time-of-flight 
between  them: 
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and  a  departing  velocity  which  is  a  function  of  current  position,  next  position,  and 
the  time-of-flight  between  them 

Vi  ~  fb  1  )  t^i+1 )  ^i+1 1  ) 

^i+1) 

The  particular  form  of  these  equations  and  method  of  evaluation  will  depend  on 
whether  the  chief  is  in  a  circular  or  elliptical  orbit. 


3.2  The  Equations  of  Motion 

The  relative  equations  of  motion  (EOM)  form  the  basis  of  the  optimization 
algorithm  presented  in  later  sections.  Specifically,  calculation  of  relative  velocity  will 
be  crucial  to  evaluating  the  cost,  which  will  be  a  function  of  AV  (since  fuel-optimality 
is  of  great  concern).  Relative  velocity  will  also  be  key  in  evaluating  our  trajectory 
constraint  (i.e.,  staying  within  the  lobe)  since  its  initial  value  along  with  a  given 
initial  relative  position  will  completely  define  a  trajectory.  We  start  with  the  general 
equations  of  relative  motion. 

Assume  the  chief  satellite  is  in  a  closed  Keplerian  orbit  and  gravity  of  the  central 
body  is  the  only  force  of  significance.  The  equations  of  relative  motion  between  a  chief 
and  deputy  satellite  are  (derivation  in  Appendix  E) 


X  —  2i)y  —  i)y  —  i> 


;.2 


X  -|- 


To 


rl{ro  +  x) 


l-\-  e  cos  V  (1  -|-  e  cos  u)  [{xo  +  xY  +  y"^  +  ^ 


y  -|-  2i'x  -\-  i)x  —  v^y 


1  - 


z  +  i?z 


(1  +  e  cos  p)  [{to  +  xY  +  +  z‘^\  ^ 

.,3 

'  =  0 


=  0 


(1  +  e  cos  p)  [(xo  +  xY  +  y'^  +  z^]  2 


0  (3a) 
(3b) 
(3c) 
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where  e  is  the  eccentricity  of  the  chief  orbit  and  is  the  instantaneous  orbital  radius, 
a  function  of  time  (or  true  anomaly): 


gsma(1  ~  e^) 
1  +  e  cos  u 


The  derivatives  of  true  anomaly  [u]  are  (derivations  in  Appendix  D) 


p 


i) 


n{l  +  ecosi^)^ 
(l-e2)l 
— 2eh^  sin  u 
1  +  e  cos  V 


(4) 


If  we  assume  that  the  deputy  is  close  to  the  chief  satellite  in  comparison  to  the  chief’s 
instantaneous  orbital  radius 


_|_  ^2  ^ 


then  the  relative  equations  reduce  to  (derivation  in  Appendix  E) 


X  —  2uy  —  i)y  — 
y  +  2i'x  +  px  —  p^y 


3  +  e  cos  p 
1  +  e  cos  p 
ecosp 
1  +  e  cos  p 


Z  +  p‘^Z 


1 

1  +  e  cos  p 


0 


=  0 
=  0 


(5a) 

(5b) 

(5c) 


It  is  possible,  and  highly  desirable  for  this  application,  to  express  time  as  a  fraction 
of  the  chief  orbit  period  as  opposed  to  an  absolute  time.  This  allows  us  to  separate 
the  relative  equations  from  a  particular  semi-major  axis  and  /r.  Referencing  Equation 
(2) 

27r  ~ 
t  =  —t 
n 
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where  t  is  time  expressed  in  fractions  of  a  chief  orbit.  The  following  conversions  are 
made  between  the  positions,  velocities,  and  accelerations  (see  Appendix  F): 

~  ~  •  27r  ~  •• 

(•)  =  (•)  (•)  =  (•)-  (•)  =  (-)^ 

n 

where  (•)  represents  x,  y,  and  2  and  tilde  over  x,  y,  and  5  indicate  that  they  are 
functions  of  chief  orbit  period.  It  is  important  to  emphasize  that  position  does  not 
scale,  thus  the  trajectories  produced  by  these  equations  are  exactly  the  same  regardless 
of  the  chief’s  semi-major  axis  or  the  value  of  y.  The  relative  equations  of  motion  are 
now  (note  that  u  =  v) 


(6a) 
(6b) 
(6c) 

where  the  derivatives  of  the  now  scaled  true  anomaly  are  dehned  as 

27r(l  -|-  e  cos  z/)^ 

(l-e2)f 

— 87r^e(l  -|-  e  cos  z/)^  sin  u 

(1  -  e2)3 

If  the  chief  is  in  a  circular  orbit  (e  =  0),  these  equations  simplify  even  further 

9  =  2tt 
9  =  0 


27r 


v  = 


p  = 


-p  = 


n 

dTT^ 

11? 


-p  = 


X  —  29y  —  9y  —  9‘^x 
y  -h  29x  +  9x  —  9‘^y 
z  +  9^z 


1  +  e  cos  9 


3  +  e  cos  9 
1  -|-  e  cos  9 
ecosz> 

1  -|-  e  cos  9 

=  0 


=  0 
=  0 
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thus 


(7a) 

(7b) 

(7c) 


X  —  47Ty  —  127r^5;  =  0 
y  +  i.'KX  =  0 
5  +  47r^5  =  0 

These  linear  time- invariant,  constant  coefficient,  differential  equations  are  the  classical 
CW  equations  [8],  albeit  in  a  less  conventional  form.  In  this  special  case,  the  Y  axis 
of  the  LVLH  frame  is  parallel  to  the  inertial  velocity  vector  of  the  chief.  The  CW 
equations  describe  two  types  of  relative  orbits  shown  in  Figure  3.  The  closed  relative 
orbit  is  constrained  to  be  a  2x1  ellipse  that  is  centered  somewhere  along  the  Y  axis. 
The  term  “2x1”  is  a  reference  to  the  fixed  ratio  of  the  semi-major  to  semi-minor 
axes  of  the  closed  relative  orbit.  The  drifting  relative  orbit  occurs  when  the  X  offset 
parameter,  a,  (Equation  (10a)  below)  is  nonzero,  a  result  of  a  difference  between  the 
deputy  and  chief’s  semi-major  axis  leading  to  a  difference  in  periods.  The  “teardrop” 
feature  size  and  shape  as  well  as  average  distance  from  the  chief  can  be  specified  [29] . 
These  differential  equations  can  be  solved  closed- form  (derivation  in  Appendix  G). 


Figure  3:  Types  of  Relative  Trajectories 
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x{i)  =  psin(27rt  +  0)  +  a 
y{t)  =  2pcos(27rf  +  9)  —  Sirai  +  b 
z{t)  =  /psin(27rt  +  0)  +  2gpcos(27rf  + 

A  useful  variant  of  the  Z  motion  is  (Equation  (97)  from  Appendix  G) 

z{t^  '^max  COs(27rt  ~\~  0) 


(8a) 

(8b) 

(8c) 


(9) 


where  the  relative  orbital  elements  (p,  a,  b,  9, 1,  q,  5max,  and  0)  are  all  functions  of  both 
the  initial  relative  position  and  velocity.  In  this  form,  the  relative  orbit  parameters 
are  given  by 


a 

b 

P 

I 

Q 

9 


^max 


-Vo  +  4:Xo 
IT 

1  A 

Vo 

7T 


ZoXo  —  47r^Zo(a  —  Xo) 


xl  +  47r2(o  —  Xo^ 

TVZoXo  +  T^Zo{a  —  Xo) 


x^  +  47r^(a  —  Xo)'^ 


(10a) 

(10b) 

(lOc) 

(lOd) 

(lOe) 

(lOf) 

(lOg) 


it>  = 


(lOh) 
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For  the  approach  presented,  the  derivatives  of  these  equations  will  also  be  necessary, 
and  are 


x(i)  =  27rpcos(27rt  +  9)  (11^) 

y{i)  =  —  dyrp  sin(27rt  +  6)  —  Svra  (Hb) 

z{i)  =  27r/pcos(27rt  +  0)  —  47rgpsin(27rt  +  0)  (He) 

or 

z(t)  =  — 27r^ma.Y  sin(27rt  +  0)  (12) 

The  equations  above  form  the  basis  of  the  controllers  developed  in  later  sections.  Note 
that  a  continuous-thrust  controller  that  keeps  the  deputy  in  a  hxed  position  relative 
to  the  chief  is  easily  synthesized  (a  linear  quadratic  regulator  for  example).  The 
research  hypothesis,  however,  is  that  a  discrete  impulsive-thrust  controller  is  more 
fuel-optimal.  To  explore  this,  these  equations  will  be  used  piecewise,  with  continuous 
position  but  discontinuous  velocities. 


3. 2. 1  Initial  &  Final  Relative  Veloeities.  Perhaps  the  most  important  infor¬ 
mation  to  extract  from  the  above  equations  of  motion  are  the  initial  and  final  relative 
velocities.  They  are  key  to  calculating  as  well  as  initializing  trajectory  propaga¬ 
tion  to  check  for  breaches  of  the  lobe  boundary.  The  initial  relative  velocity  under 
the  circular  chief  orbit  assumption  is  (derivation  in  Appendix  I.l) 


Xo 

-45+67rf(7 

2-2(7 

0 

4S-6iTf 

-2+2(7 

0 

8-67rT5-8C 

8-67rT5-8(7 

8-6nfS-8C 

8-6nfS-8C 

Vo 

=  27r 

-14+127rf5+14C 

-S 

0 

2-2(7 

S 

0 

8-67rT5-8C 

8-6nfS-8C 

8-67rT5-8(7 

8-6nfS-8C 

Zo 

0 

0 

c 

s 

0 

0 

1 
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Mo 


Xo 

Vo 

Xf 

Vf 

% 


(13) 
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where 


T 

^  =  p 


S  =  sin(27rT')  C  =  cos(27rT') 


and  Mo  is  the  transformation  matrix  that  maps  the  initial  and  final  positions  into 
initial  relative  velocity.  The  final  relative  velocity  is  (derivation  in  Appendix  1.2) 


Xf 

-4S'+67rT 

-2+2C 

0 

4S-67tTC 

2-2C 

8-6iTfS-8C 

8-67rfS-SC 

8-67rf5-8C 

8-67rfS-SC 

Vf 

=  27r 

2-2C 

-S 

0 

-14+127rf;§+14C 

S 

S-GwfS-SC 

8-6nfS-8C 

8-67rfS'-8C 

8-6nfS-8C 

0 

0 

1 

s 

0 

0 

Mf 


Xo 

Vo 

Vf 


(14) 

where  Mf  maps  initial  and  final  positions  into  final  relative  velocity.  Note  that  there 
are  two  singularities  in  Mq  and  M/  at 


sin(27rT')  =  0 

8  —  GyrT  sin(27rT')  —  8cos(27rT')  =  0 


The  first,  sin(27rT),  is  easily  solved  for  time-of-fiight:  T  =  0,  -  and  is  only 

active  if  there  is  Z  motion.  Physically  the  deputy  is  passing  through  the  relative 
orbit’s  axis  of  rotation.  This  is  the  intersection  of  an  infinite  number  of  orbits  and 
thus  the  linear  system  is  indeterminate. 

The  second  singularity  (see  Figure  4)  is  not  as  easy  to  find  except  for  the  first  two 
zeros:  T  =  0  and  T  =  1.  If  not  for  the  secular  term,  a  closed- form  expression  could  be 
found.  Unfortunately  we  must  find  further  zeros  numerically.  Figure  5  illustrates  the 
change  in  the  interval  between  zero  crossings,  an  interval  that  is  clearly  converging  to 
0.5  chief  orbits.  Note  that  the  zero  crossing  number  is  an  integer  value,  however,  the 
graph  points  are  connected  for  clarity.  There  appears  to  be  two  exponential  decay 
patterns:  one  for  even  and  one  for  odd  crossings.  It  may  not  be  difficult  to  express 
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these  decays  closed-form,  however,  we  are  normally  only  interested  in  times-of-flight 
between  0  and  1,  therefore  finding  the  zero  crossing  pattern  is  not  germane  to  this 
research. 


Figure  4:  Plot  of  f{T)  =  8  —  QwTS  —  8C 


Figure  5:  Time  Between  the  Zeros  of  f{T)  =  8  —  GttTS  —  8C 


3.2.2  Equilibrium  Points  of  the  EOM.  It  is  always  prudent  to  find  the 
equilibrium  points  of  any  set  of  differential  equations.  This  is  especially  important 
for  this  analysis  since  equilibrium  represents  zero-fuel  solutions  and  may  be  the  ideal 
place  to  put  the  deputy  satellite.  While  in  equilibrium,  fuel  is  used  only  to  reject 
disturbances  and  linearization  errors  in  order  to  stay  at  the  equilibrium  point.  Starting 
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with  the  CW  equations  of  motion  (Equation  (7)),  let 


x=y=z=x=y=z=0 


then 


X  —  A'Ky  —  12t^‘^x  (0)  —  47r(0)  —  127r^5;  =  — 127r^5;  =  0 

y  +  Attx  (0)  +  47r(0)  =  0 

5  +  Att'^z  (0)  +  Atv^z  =  0 


Thus  X  =  0  and  z  =  0  without  any  restriction  on  y,  meaning  the  entire  Y  axis  is  an 
equilibrium  point  (under  linear  assumptions).  Similarly,  the  equilibrium  condition  for 
the  non- linearized  equations  of  motion  can  also  be  found.  Applying  the  zero  derivative 
conditions  to  Equation  (74),  Appendix  E. 


(0)  -  (0)  -  n2(ro  x) 
(0)  -S  (0)  -  n'^y 
(0)  -|-  n^z 


1  - 


[(ro+x)^+y^+z^]  2 
^3 


[{ro+x)‘2W+^‘^Y 

^3 


[(ro+x)2+2/2+^2]^ 


=  0 


=  0 
0 


Note  there  is  no  longer  any  advantage  to  using  the  scaled  tilde  versions  of  the  relative 
position  and  velocity  since  To  appears  in  these  equations.  Since  To  and  n  are  always 
positive  and  (vo  +  xY  +  y'^  Y  z'^  non-zero,  it  is  clear  that  z  =  0.  In  the  other  two 
directions,  the  term 

^3 

1 - 2 - - 

[{vo  +  xy  +  y‘^  +  z'^y 
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going  to  zero  will  satisfy  the  equilibrium  conditions  (so  will  x  =  y  =  0,  the  trivial 
solution).  Solving  yields 


3 

[vo  +  +  y‘^  +  z'^]^  -  rl  =  0 

[{ro  +  X)"^  +  _|_  ^2j2  =  ^3 

[{to  +  xY  +  +  z^]  =  rl 

Note  that  the  left  hand  side  is  the  square  of  the  magnitude  of  the  deputy’s  inertial 
position  (M)  written  in  relative  coordinates.  Thus 

\\M  =  rl  ^  ||M||2  =  r, 

Therefore  the  deputy  must  be  in  same  size  circular  orbit  as  the  chief.  Also,  since 
z  =  0,  the  deputy’s  orbit  must  be  co-planar  with  the  chief  satellite.  All  together,  this 
means  the  chief  and  deputy  share  the  same  orbit  but  may  have  different  arguments 
of  latitude. 

3.2.3  Symmetries  of  the  Circular  Chief  Relative  EOM.  The  circular  chief 
equations  of  motion  contain  symmetries  that  allow  us  to  generalize  results  from  lobes 
in  one  quadrant  of  the  relative  frame  to  results  in  any  of  the  other  four  quadrants, 
thereby  signihcantly  reducing  the  number  of  results  we  need  to  produce  and  examine. 
First  recall  that  a  function  is  odd  if  [24] 

f(-x)  =  -f{x) 
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The  equation  for  initial  velocity  and  final  velocity  are  obviously  odd  since  they  are 
linear  in  the  states: 
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Converting  Equation  (108),  Appendix  I  to  time  units  of  chief  orbit  period: 
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where  S  =  sin(27rt)  and  C  =  cos(27rt).  As  a  function  solely  of  starting  and  ending 
relative  positions  and  time-of-fiight: 
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Once  again  we  have  a  function  that  is  linear  in  the  states  and  thus  an  odd  function. 
An  example  is  provided  in  Figure  6;  the  deputy  starts  at  (2  km,  1  km)  and  (-2  km, 
-1  km)  and  ends  at  (1  km,  0.25  km)  and  (-1  km,  -0.25  km)  both  with  a  time-of-flight 
of  0.35  chief  orbits.  Each  point  on  one  trajectory  can  be  rotated  about  the  origin  tt 
radians  while  maintaining  equidistance  from  the  origin  to  get  the  corresponding  point 
on  the  other  trajectory.  These  two  points  occur  at  the  same  i  in  their  respective 
trajectories. 


Figure  6:  Symmetry  of  the  Circular  Chief  Equations  of  Motion 


3.2.4  Error  Analysis.  Several  different  forms  of  the  equations  of  motion  have 
been  presented.  Some  have  linearizing  assumptions  and  all  do  not  model  perturbations 
and  other  error  sources.  Before  posing  and  solving  the  optimization  method,  the 
validity  and  accuracy  of  these  equations  must  be  verihed. 

As  with  most  engineering  problems,  the  optimization  algorithm  is  based  on  a 
reduced-order  model  for  which  mathematically  tractable  solutions  are  available.  As 
we  will  see,  the  short  total  time-of-flight  associated  with  this  problem  keep  total  errors 
within  a  reasonable  range.  There  are  four  main  sources  of  error: 

1.  Integration  errors 
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2.  Numerical  errors  during  the  conversion  between  inertial  and  relative  reference 
frames 

3.  Linearization  error  from  the  closeness  approximation 

4.  Unmodeled  perturbations 

Integration  errors  during  trajectory  propagation  are  a  standard  problem  with  orbital 
research  and  can  be  mitigated  with  the  proper  choice  of  a  maximum  propagation 
timestep.  Figure  7  displays  the  hnal  position  error  for  trajectories  produced  with 
Equation  (7)  (CW  solutions  as  functions  of  chief  orbit  period)  for  various  maximum 
timesteps.  The  error  is  calculated  by  comparing  the  hnal  position  of  each  integration 
run  over  one  chief  period  versus  an  extremely  small  timestep  run  (1  x  10“^).  In  order 
to  achieve  integration  errors  less  than  1  x  10“®  meters,  0.003  has  been  chosen  as  the 
largest  acceptable  timestep. 


Figure  7:  Trajectory  Integration  Errors  Over  One  Chief  Period 

Additional  numerical  errors  are  encountered  when  using  the  inertially  propa¬ 
gated  truth  model  which  requires  the  rotation  of  inertial  relative  position  and  velocity 
vectors  into  the  LVLH  frame  (Appendix  H).  These  results  are  differenced  from  trajec¬ 
tory  propagations  of  Equation  (3)  with  e  =  0  (which  yields  u  =  n  and  h  =  0).  Since 
there  are  no  linearization  assumptions,  errors  between  the  trajectories  are  numeri- 
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cal  in  nature  caused  by  an  additional  mathematical  operation  within  the  algorithm 
during  which  roundoff  and  truncation  errors  occur.  There  is  little  that  can  be  done 
to  mitigate  this  error  besides  proper  maximum  time  step  choice  as  discussed  above. 
Figure  8  shows  the  magnitude  and  trends  of  the  transformation  errors.  Although  it 


Figure  8:  Numerical  Error  Due  to  Inertial  to  Relative  Frame  Rotations 

is  clear  that  the  numerical  imprecision  builds  with  time,  at  least  in  the  Y  direction, 
the  errors  are  well  within  acceptable  limits  even  after  one  chief  orbit.  Although  errors 
are  not  scaled  to  a  p  to  To  ratio.  Figure  8  shows  the  results  of  four  separate  runs 
with  different  values  of  p.  To  summarize,  the  first  two  error  sources  are  standard 
assumptions  for  orbital  mechanics  problems.  Most  missions  to  which  this  technique 
will  be  applied  are  of  short  duration  and  close  proximity,  both  of  which  work  in  our 
favor  to  reduce  errors. 

The  third  source  of  error  occurs  when  we  make  the  following  linearization  as¬ 
sumption  for  model  simplification  (see  Appendix  E): 

x‘^ 

—  fs  —  —  PS  0  and  To  -|-  3a:  ~ 


where  is  the  instantaneous  orbital  radius  of  the  chief.  We  naturally  expect  to  see 
errors  increase  as  average  chief  to  deputy  distance  increases  and/or  as  the  orbital 
radius  decreases.  The  relative  orbit  element  p  is  a  good  parameter  with  which  to 
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measure  average  distance  from  the  chief.  Figure  9  shows  errors  for  a  deputy  in  a 
closed  relative  orbit  about  the  chief  (i.e.  the  X  offset,  a,  is  zero).  Errors  for  closed 


Chief’s  Semi-Major  Axis  =  6628.14  km 


Figure  9:  Linearization  Error  in  a  Closed  Relative  Orbit  (a  =  0  km) 


relative  orbits  are  nearly  periodic  in  the  X  axis  but  have  a  strong  secular  drift  in  the 
Y  axis.  Although  setting  a  =  0  cancels  the  first-order  secular  term,  it  leaves  higher- 
order  secular  terms  that  build  in  error  with  time.  Eventually,  the  deputy  will  drift 
far  enough  away  from  the  chief  to  invalidate  the  linearized  equations  of  motion.  This 
happens  even  more  quickly  in  the  drifting  relative  orbit  case  (Figure  10)  in  which 
the  substantial  first-order  drift  term  adds  to  the  growth  in  distance  between  chief 
and  deputy.  It  is  clear  that  a  short  time-of-fiight  for  each  trajectory  leg  is  key  to 
the  validity  of  the  reduced-order  model  upon  which  the  optimization  routine  is  built. 
Figure  11  shows  the  maximum  error  of  the  deputy  after  one  period  of  the  chief  orbit 
due  to  linearization  errors. 

The  final  source  of  error  is  due  to  unmodeled  perturbations.  The  two  largest 
perturbations  are  due  to  the  oblate  Earth  (J2)  and  atmospheric  drag  [59].  Both 
are  inversely  proportional  to  semi-major  axis  and  are  thus  stronger  at  lower  alti¬ 
tudes.  At  these  lower  altitudes,  the  Earth  is  better  modeled  as  non- spherical  and 
with  non-uniform  density,  leading  to  a  infinite  series  of  correction  terms  to  the  two 
body  problem.  The  first  and  most  significant  of  these  terms  is  J2  (derivation  and 
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Chief’s  Semi-Major  Axis  =  6628.14  km  Chief’s  Semi-Major  Axis  =  6628.14  km 


Figure  10:  Linearization  Error  in  a  Drifting  Relative  Orbit  (a  =  2  km) 


Orbital  Altitude  (km)  p  (km)  Orbital  Altitude  (km)  p  (km) 


(a)  Closed  Relative  Orbit  (a  =  0  km)  (b)  Drifting  Relative  Orbit  (a  =  2  km) 

Figure  11:  Maximum  Linearization  Errors  Over  One  Chief  Period 


discussion  in  Appendix  0.2): 
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where  is  the  radius  of  the  Earth,  X,  Y,  and  Z  are  the  inertial  coordinates  the 
satellite  and 


d  =  VV2  +  y2  ^ 
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Measurements  of  the  zonal,  tesseral,  and  sectorial  coefficients  reveal  that  J2  is  at  least 
400  times  larger  than  the  next  most  signihcant  term  and  is  thus  the  only  harmonic 
considered  here.  Figure  12  shows  the  error  between  the  full  nonlinear  CW  model  and 
a  truth  model  including  J2  effects  at  an  orbital  altitude  of  250  km. 


Chief’s  Semi-Major  Axis  =  6628.14  km 


Time  (Chief  Orbit  Fractions) 


Chief’s  Semi-Major  Axis  =  6628.14  km 


Figure  12:  Perturbations  Due  to  J2  (Orbital  Altitude  =  250  km) 


The  second  source  of  error  considered  is  atmospheric  drag,  modeled  as  a  force 
opposing  the  relative  wind  of  the  satellite  and  based  on  an  exponential  model  of 
atmospheric  density  (derivation  and  discussion  in  Appendix  0.3). 


^  K-el 

®drag  TTd  Po®  ^  ^rel 
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(16) 


where  V^ei  is  the  satellite’s  velocity  relative  to  the  atmosphere,  Be  is  the  ballistic 
coefficient,  po  is  the  nominal  density,  ho  is  the  reference  altitude,  hd  is  the  satellite’s 
altitude,  and  H  is  the  scale  height.  The  error  charts  presented  in  Figure  13  are  a  worse 
case  scenario  of  low  altitude  (250  km)  and  disparate  ballistic  coefficients  between  chief 
and  deputy.  Representative  values  of  the  ballistic  coefficient  are  chosen  from  historical 
satellites,  a  micro  satellite  for  the  depnty  and  a  larger  scientihe  class  satellite  for  the 
chief. 

Chief  Be  =  25  ^  Deputy  Be  =  128  ^ 
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Chief’s  Semi-Major  Axis  =  6628.1 4  km  Chief’s  Semi-Major  Axis  =  6628.1 4  km 


Figure  13:  Perturbations  Due  to  Atmospheric  Drag  (Orbital  Altitude  =  250  km) 


Assuming  proper  choice  of  time  step,  the  largest  sources  of  error  by  far  are 
linearization  and  J2.  All  of  the  results  presented  in  this  research  are  for  lobes  that 
are  less  than  3  km  from  the  chief;  therefore,  over  one  chief  orbit  period  the  total 
error  should  be  less  than  100  meters.  Of  course,  most  hovering  trajectories  will  have 
a  time-of-flight  well  under  a  chief  orbit  period,  therefore  only  minor  changes  to  the 
path  will  be  noticeable. 


3.3  The  Cost  Function 


The  problem  statement  motivates  an  optimization  problem  that  requires  the 
formulation  of  a  cost  function.  Ideally,  we  want  to  maximize  the  time  spent  inside 
the  lobe  per  unit  of  fuel  expended,  or  equivalently  we  can  minimize  the  cost  function 
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where  AVj  is  the  instantaneous  change  of  velocity  required  at  the  ith  point  of  impulsive 
thrust,  k  is  the  number  of  legs  (a  leg  being  defined  as  the  trajectory  between  two 
impulsive  thrusts),  and  the  subscript  F  (for  final)  represents  a  possible  exit  burn 
and  time-of-flight.  AV  is  calculated  by  taking  the  Euclidean  norm  of  the  difference 


31 


between  pre-  and  post-thrust  relative  velocities  (derivation  in  Appendix  K). 


A14  =  II  Ai;i||2  =  \\vl  -  Vi  II2  =  ^ {it  -  )^  +  {vt  -  %  Y 

Using  AU^  yields  a  more  mathematically  compact  solution  for  impulsive  thrust  and 
allows  for  closed- form  solutions  of  the  cost  functions’s  gradient.  This  gradient  yielded 
only  minimal  performance  gains  during  calculation  of  the  optimal  trajectory  and 
therefore  was  not  implemented  in  the  final  results;  however,  the  derivation  of  the 
gradient  is  provided  to  future  researchers  in  Appendix  P.  Since  AU  is  always  positive, 
minimizing  the  square  will  minimize  the  value  itself,  thus  our  cost  function  is 

j  _  +  AU^^  +  IWj  +  ■■■  +  +  AU| 

^1,2  +  2^2,3  +  ^3,4  +  •••  +  Tk,k+l  +  Tp 

Impulsive  maneuvers  and  the  time-of-flight  can  be  expressed  in  a  wide  range  of 
units.  Initial  investigations  into  the  problem  found  unacceptably  flat  cost  functions, 
thus  encouraging  the  proper  use  of  scaling.  Experimentation  with  this  cost  function 
showed  that  it  was  generally  desirable  to  keep  both  the  sum  of  the  AU’s  and  the  total 
time-of-flight  somewhere  in  the  range  of  0  to  1.  For  the  numerator,  an  appropriate 
scaling  term  is  the  continuous  AU  (AVc)  required  to  keep  the  deputy  at  the  x^in 
position,  which  is  the  optimal  place  to  stay  for  a  continuous  thrust  scenario.  As  will 
be  shown  in  Section  4.6,  the  minimum  continuous  thrust  AU  is 

AUc  =  Gx^inmvTT 

where  Tp  is  the  total  time-of-flight.  Likewise,  time-of-flight  is  conveniently  scaled  by 
the  chief’s  orbit  period  (P): 

p  [^1,2  +  ^2,3  +  ^3,4  +  •••  +  T'fc,fc+1  +  Tp]  =  Ti^2  +  ^2,3  +  ^3,4  +  •••  +  T/c,/c+l  +  Tp 
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Combining  the  scaling  of  AV  and  time-of-flight,  the  cost  function  is  now 


^  [AV^^  +  A\q  +  AVi  +  ...  +  AVi  +  AV}] 

Ti,2  +  ^2,3  +  T^,A  +  •••  +  Tk^k+1  +  Tp 


As  discussed  in  Section  3.2,  it  is  desirable  to  have  equations  that  are  not  a  function 
of  the  chief’s  semi-major  axis  or  /r.  To  that  end,  define  specific  delta  V,  AV  as  (see 
Appendix  K.2) 


AC  = 


AC  AC 


(17) 


n  2w 

where  AC  is  the  Euclidean  norm  of  the  difference  of  the  relative  velocities  expressed 
with  time  in  units  of  chief  orbit  period.  The  equation  for  the  square  of  AC  is  (deriva¬ 
tion  in  Appendix  K) 
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with  R  defined  as 
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where 
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T-  =  r+  = 

S~  =  sin(27rf'j_i_j)  5+  =  sin(27rT;^j+i) 

c~  =  cos(27rf;_i,j)  (7+  =  cos(27rTj,i+i) 

and  the  positions  are  (derivation  in  Appendix  B) 


Xi  =  j  cos  a  sin  (3  + 
Vi  =  J  sin  a  sin  f3  + 


T^Ty  cos 


^ cos2  {ipi  -ri)+T^  sin^  (^j  -  7]) 

_ TxTy  sin  xpi _ 

cos2  {i/ji  -r])  +  T^  sin^  {ipi  -  rj) 


Noting  that  the  continuous  thrust  scaling  term  as  a  specific  AV  is 


=  ^  =  36xi,y-n 


This  results  in  a  cost  function  for  the  problem  posed  herein  of 
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3.4  Interior  Thrust  Points 

As  a  simplifying  assumption,  all  thrusting  is  required  to  occur  on  the  lobe 
boundary.  This  removes  a  degree  of  freedom  from  the  problem  and  thus  only  a 
single  coordinate  (ip)  is  required  to  specify  a  thrust  point.  In  addition  to  being  a 
natural  assumption  to  start  with,  it  turns  out  that  the  cost  function  is  ill-posed  to 
handle  interior  points.  Let’s  assume  a  minimum-cost  single  leg  trajectory  is  found. 
Without  additional  constraints  on  the  lower  bound  of  the  time-of-flight  for  each  leg 
or  a  minimum  AN  magnitude,  there  is  nothing  to  prevent  the  algorithm  from  simply 
parsing  that  single  leg  solution  into  as  many  “legs”  as  are  requested.  Application 
of  those  additional  constraints  or  a  reposed  cost  function  that  is  more  favorable  to 
interior  points  is  left  to  future  researchers. 
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3. 5  Units 


The  non-dimensional  time,  T,  is  of  course  unitless  and  more  importantly,  not 
tied  to  a  particular  orbital  radius.  A  similar  disentanglement  between  /SV  and  orbital 
radius  is  made  possible  with  specihc  /SV  and  has  the  following  units 

ATT  length 

=  length 

n  tv-A 

time 

Likewise  the  units  of  relative  velocity  when  time  is  expressed  in  units  of  chief  orbital 
period  are 

^  .  length 

(•)  =  -(-)  =  ^  =  length 

time 

and  the  relative  accelerations 

4^2  ..  Isilgth 

(■)  =  —(-)  =  -^  =  length 

Since  ^  are  scaled  by  AVq,  the  numerator  of  J  is  unitless.  Dividing  by  the 
unitless  total  time-of-flight  we  note  that 

J  =  unitless 


as  desired. 

3. 6  The  Constraint 

From  a  mission  planner’s  standpoint,  the  only  constraint  on  the  deputy  satel¬ 
lite’s  motion  is  that  it  stay  within  the  prescribed  lobe.  While  there  are  several  methods 
to  pose  this  constraint,  one  relatively  simple  way  is  to  force  the  time-of-flight  to  be 
smaller  than  or  equal  to  some  maximum  time-of-flight.  This  maximum  time-of-flight 
is  naturally  defined  as  as  the  largest  time-of-flight  for  which  the  deputy’s  entire  tra¬ 
jectory  remains  inside  the  lobe.  If  two  positions  are  chosen  in  the  relative  frame  along 
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Figure  14:  Notional  Time-of-Flight  Comparisons 

with  a  time-of-flight  between  them,  the  relative  velocities  can  be  found  via  Equations 
(13)  and  (14).  An  appropriately  small  time-of-flight  produces  a  trajectory  that  closely 
approximates  a  straight  line  between  the  two  points  and  has  a  large  magnitude  rel¬ 
ative  velocity.  As  the  time-of-flight  is  increased,  the  trajectory  exhibits  larger  and 
larger  loops  that  may  or  may  not  cross  back  over  themselves  (teardrop  maneuvers). 
As  demonstrated  in  Figure  14,  there  is  a  maximum  time-of-flight  after  which  any 
larger  time-of-flight  creates  a  trajectory  that  partially  leaves  the  constraint  lobe.  Any 
time-of-flight  smaller  than  this  maximum  will  satisfy  the  mission  planner’s  constraint. 
Thus  the  constraint  is  formulated  as 

(f/'l ,  -02)  (21) 

The  maximum  time-of-flight  between  any  two  points  is  precomputable,  and  a  surface 
is  easily  generated  for  a  speciflc  lobe  size  and  shape  (example  Figure  15).  Interpolating 
between  points  is  a  very  effective  way  to  calculate  the  constraint  during  optimization 
searches.  Note  that  unless  otherwise  stated,  all  constraint  surfaces  in  this  document 
use  the  same  normalized  colorbar.  Thus  green  in  one  flgure  is  the  same  value  of  T^ax 
as  green  in  another  figure. 

All  constraint  surfaces  have  a  valley  of  Tmax  =  0,  the  set  of  -^’s  from  which 
the  deputy  cannot  start  and  end  without  leaving  the  lobe  no  matter  how  small  the 
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Figure  15:  Maximum  Time-of-Flight  {a  =  45°,  7  =  2  km,  Tj,  =  1  km,  =  0.5  km) 

time-of-flight.  These  points  are  located  on  the  side  of  the  lobe  closest  to  the  Y  axis 
(see  Figure  16).  Since  our  lobe  is  a  closed  concave  shape,  we  can  always  draw  a 
straight  line  between  any  two  points  that  stays  completely  within  the  lobe  (and  as 
noted  previously,  as  T^ax  0  the  trajectory  becomes  a  straight  line),  therefore  this 
valley  will  always  be  along  the  line  'tpi  =  ip2-  Lobes  that  do  not  intersect  the  Y  axis 
will  also  have  a  peak  formation.  Note  that  this  surface  repeats  in  both  the  and  xp2 
directions  since  the  angle  is  modulo  27r  and  thus  the  four  peaks  in  Figure  15  are  in 
fact  a  single  peak. 

3.6.1  Calculation  ofT^ax-  Unfortunately,  no  closed- form  solution  for  Tjnax  is 
available,  prompting  a  numerical  solution.  Two  methods  of  calculation  were  explored 
for  this  research.  The  most  robust  and  most  computationally  expensive  method  is  to 
increment  time-of-flight  and  check  the  resulting  trajectory,  converging  on  a  solution 
when  the  time-of-flight  produces  a  trajectory  that  leaves  the  lobe.  In  practice,  the 
midpoint  of  an  upper  and  lower  bound  is  used  for  the  trajectory  check.  If  the  trajec¬ 
tory  stays  within  the  lobe,  the  midpoint  is  now  the  new  lower  bound  while  it  becomes 
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Figure  16:  Points  of  Zero  Maximum  Time-of-Flight 

the  new  upper  bound  if  it  violates  the  boundary.  The  flowchart  of  this  process  is 
found  in  Figure  17. 

The  second  numerical  technique  involves  identifying  two  conditions  of  a  trajec¬ 
tory  which  has  a  single  point  of  tangency  with  the  lobe.  The  first  is  that  the  point  of 
tangency  (^t)  occurs  on  the  lobe  boundary: 

rii^T)  -  !  =  0 

J cos2  {^pT  -  h)  +  F?  i'ipT  -  v) 

and  the  second  that  the  relative  velocity  of  trajectory  at  'ipT  is  tangent  to  the  lobe: 

xt  cos  ipT  +  Vt  sin  ipr  =  ^ 

These  two  nonlinear  equations  must  be  solved  numerically,  but  algorithms  such  as 
FSOLVE  in  Matlab®  do  the  job  quickly  and  efficiently.  This  method  was  found  to 
be  less  robust,  however,  when  the  point  of  tangency  occurred  at  or  near  the  beginning 
or  end  point  of  the  trajectory.  Because  accuracy  was  valued  over  speed  for  this 
particular  application,  the  first  method  was  chosen  to  produce  the  constraint  surface. 
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Figure  17:  Flowchart  for  Finding  Tmax 

Definition  and  calculation  method  in  hand,  we  turn  next  to  several  properties 
of  the  constraint  surface  that  will  narrow  the  relevant  lobe  parameter  space  and  yield 
more  efficient  calculations  of  the  constraint  surface. 

3.6.2  Invariance  of  the  Constraint  Surface  to  pl-  Under  the  CW  assump¬ 
tions,  it  can  be  shown  that  the  constraint  surface  is  invariant  to  the  y  coordinate  of 
the  lobe  center  {yh)-  The  first  indication  of  this  is  found  in  the  original  differential 
equations  (Equation  (7)),  none  of  which  are  functions  of  the  y  coordinate.  However, 
clearer  proof  is  found  by  examining  the  initial  and  final  relative  velocity  equations 
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(Equations  (13)  and  (14)).  Let  D  =  8  —  GttT sin(27rr)  —  8cos(27rr),  then: 


Likewise, 


Xo  = 


Vo 


,  ,  ,  2-2  cos(27rT)  -2  +  2  cos(27rr) 

fo[Xo,  Xf)  +  - - - Vo  +  - ^ - Vf 

/.(*.,*/)  -  -  y.) 

^  sin(27rT')  ^  sin(27rT')  ^ 

QoyXoiXf)  —  yo-\  ^  Vf 


D 


D 


f~  -sin(27rT)  ^ 

=  go{Xo,  Xf) - - - [yf  -  yo) 


=  ffiS^cXf)  -  yo) 

~  /-  ~  \  sin(27rf’)  _ 

Vf  ~  9f{xo^Xf)  —  {yf~yo) 


We  see  that  the  relative  velocities  are  functions  of  Ay  only  and  not  the  absolute 
positions  causing  the  constraint  surface  to  remain  unchanged  as  it  slides  up  or  down 
a  line  parallel  to  the  Y  axis  (see  Figure  18).  This  is  convenient  for  users  desiring  to 
precompute  the  constraint. 
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o 

Invariance  of  the  Constraint  Surface  to  y^ 


3.6.3  Symmetries  Between  Constraint  Surfaees  in  the  Right  and  Left  Hand 
Planes.  In  Section  3.2.3,  the  equations  of  motion  were  shown  to  be  odd  functions 
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meaning  that  every  trajectory  has  a  twin  trajectory  in  the  opposite  quadrant  that 
is  rotated  by  vr  (as,  for  example,  Figure  19).  Taking  this  concept  a  step  further,  if 
we  project  our  lobe  into  the  opposite  quadrant  as  well,  we  hnd  that  the  maximum 
time-of-flight  (Tmax)  of  a  trajectory  between  and  'ip2  in  the  1st  quadrant  is  also  the 
maximum  time-of-flight  between  'tpi  +  %  and  -^2  +  tt  in  the  3rd  quadrant. 


Tmax(t/^l)  1/^2)  Tkiax(t^l  T  1^2  T 


Lobe  1 


Lobe  2 


It  is  the  symmetry  of  the  lobe  that  allows  this  transformation  to  occur  (i.e.,  our 


-3000  -2000  -1000  0  1000  2000  3000 

X  Axis  (m)  -  Radial  Direction 


Figure  19:  Symmetric  Lobes  and  Reflective  Trajectories 


mirror  ellipse  orientation  is  77  =  rj  +  n).  The  two  trajectories  are  negatives  of  each 

other: 


x{i)  =  —x(i) 
y{i)-  =  -y{t) 
for  0  <  i  <  Tmax 
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To  prove  this,  observe  that  the  second  lobe  is  located  at  a  =  «  +  tt.  Then 


X  =  7COs(q;  +  tt)  sin/3  + 
y~  =  7sin(Q;  +  tt)  sin/3  + 


T^Ty  COs(^j  +  Tt) 


Noting  that 


^ cos2  -  [?7  +  tt])  +  sin^  -  [?7  +  tt]) 
_ TxTy  +  tt) _ 

y r2  cos^  -  [r/  +  tt])  +  rj  sin^  (?/^i  -  [^  +  t’’]) 


sin('^j  +  tt) 

=  —  sin  xpi 

cos{ipi  +  tt) 

=  —  cos  Ipi 

sin^(^i  -[r) +  7:]) 

=  -  T]) 

cos^ {'ipi  -  [r/  +  tt]) 

=  cos^  {ijji  -  y) 

then 


a:  = 


y  =  - 


7  cos  a  sin  /3  — 
7  sin  a  sin  /3  — 


ra;r„  cos 


y r2  cos2  -r])  +  T^  sin^  {ijji  -  rj) 

_ TxTy  sin  -ipi _ 

y cos2  (?/;j  -  r^)  +  sin^  -  77) 


=  —a: 


=  -y 


Thus  we  see  that  transforming  the  constraint  surface  from  the  1st  to  the  3rd  quadrant 
involves  a  shift  of  the  starting  and  ending  ?//  by  tt  radians.  Applying  the  invariance 
to  the  yi^  discussed  in  Section  3.6.2,  we  can  translate  the  negative  lobe  into  the  2nd 
quadrant  (illustrated  in  Figure  20).  So,  for  the  computational  price  of  one  lobe  we 
get  a  family  of  lobes  along  lines  parallel  to  the  Y  axis  in  both  the  right  and  left  hand 
planes  (Figure  21).  Lobes  are  therefore  functions  only  of  the  x  coordinate  of  their 
centers  (x^),  size  {Tx,Ty),  and  orientation  (rj).  The  AV  calculations  are  also  functions 
of  Ay  as  they  are  derived  directly  from  the  relative  velocity  equations.  Thus,  with 
the  constraint  surfaces  and  AV  calculations  being  invariant  to  being  placed  in  either 
the  right  or  left  hand  planes  (same  distance  from  the  Y  axis)  and  identical  in  the 
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Figure  20:  Equivalent  Constraint  Surfaces  (a  =  135°, 7  =  2  km,Ta;  =  1  km, = 
0.5  km) 


upper  and  lower  half  planes,  we  can  conclude  that  optimal  trajectories  found  in  lobes 
in  the  1st  quadrant  are  indicative  of  lobes  in  any  other  quadrant.  We  will  therefore 
concentrate  solely  on  1st  quadrant  lobes  in  the  results  chapter,  as  these  results  can 
easily  be  mapped  to  any  quadrant. 


3.6.4  Lobes  Symmetric  About  Their  Horizontal  Axis.  Lobes  that  are  sym¬ 
metric  about  their  horizontal  axes  have  an  additional  symmetry  in  their  constraint 
surface  that  can  halve  the  needed  computation  time.  For  these  lobes, 


rmax(t/'l,  271  -  ^2)  =  T^a.xi'4’2,  StT  -  ^i) 


An  example  is  provided  in  Figure  22. 

3.6.5  Lobes  that  Intersect  the  Y  Axis.  Section  3.2.2  demonstrated  that  the 
Y  axis  is  the  set  of  equilibrium  points  for  the  linearized  equations  of  motion  (chief 
in  a  circular  orbit).  We  expect  that  lobes  containing  the  Y  axis  will  have  unique 
properties.  An  infinite  number  of  2x1  closed  relative  orbit  ellipses  can  be  placed  within 
lobes  that  have  two  points  of  intersection  with  the  F,  and  certain  combinations  of  rj 
and  xl  allow  a  single  closed  relative  orbit  to  be  tangent  to  both  points  of  intersection. 
This  provides  an  opportunity  to  start  and  end  at  either  the  upper  or  lower  Y  axis 
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Figure  21:  Equivalent  Lobes 

intersection  point  and  come  back  to  the  same  point  in  one  full  chief  orbit.  It  also 
provides  a  trajectory  that  in  the  absence  of  perturbations  and  ignoring  linearization 
error  can  be  maintained  with  zero  fuel.  Note  that  the  relative  velocity  (Equations 
(13)  and  (14))  is  undefined  at  T  =  1,  a  singularity  of  which  we  must  be  cognizant 
when  evaluating  these  orbits.  If  the  upper  lobe  intersection  point  is  'ifju  and  the  lower 


intersection  xpL,  then  there  are  four  regions  about  which  there  is 

trajectory  that  is  or  nearly  is  one  full  chief  orbit;  they  are: 

an  opportunity  for  a 

(1) 

(22a) 

(2) 

(22b) 

(3) 

(22c) 

(4) 

(22d) 
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Figure  22:  Lobe  Symmetric  About  Horizontal  Axis  {a  =  0°,7  =  1.414  km,  r^;  = 

1  km,  Ty  =  0.5  km,  rj  =  0°) 


Pairs  (1)  and  (4)  yield  2x1  closed  relative  orbit  ellipses  (Figure  23(a))  that  have  a 
maximum  possible  semi-major  axis  of  (derivation  in  Appendix  M.l) 


P  < 


r{'i/^u)  siii'ipu  —  r{'i/^L)  siii'ipi 


(23) 


Pairs  (2)  and  (3)  represent  drifting  relative  orbits  with  small  p  and  a  values  that 
are  nearly  contained  within  the  lobe  over  one  chief  orbit  period  (see  Figure  23(b)). 
The  constraint  surface  for  lobes  that  do  not  intersect  the  Y  axis  have  a  single  peak 
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Figure  23:  Trajectories  of  Nearly  One  Chief  Orbit  Period 


formation. 


As  the  lobe  passes  through  the  Y  axis,  this  single  peak  splits  into  four 
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peaks  attached  to  the  four  conditions  in  Equation  (22).  Peaks  surrounding  angle 
pairs  (1)  and  (4)  will  slide  along  the  line  'ipi  =  xp2  while  (2)  and  (3)  will  slide  along 
2tt  —  ipi  =  '^2,  all  corresponding  to  the  ip^s  of  lobe’s  intersection  with  the  Y  axis. 
To  illustrate  this,  Figure  24  displays  the  constraint  surface  and  peak  location  for  the 
following  lobe: 


a  =  90°  =  2  km  ry  =  1  km  Ty  =  0.5  km  rj  =  20°  (24) 

thus  the  center  is  located  qXxl  =  Q  and  the  intersection  points  are  in  the  neighborhood 
of 
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Figure  24:  Lobe  Intersection  of  the  Y  Axis  (a  =  90°,  7  =  2  km,  Ta,  =  1  km,  = 
0.5  km,  r]  =  20°) 


3. 6. 6  Sensitivity  of  the  Constraint  Surfaee  to  Lobe  Size.  Larger  lobes  provide 
more  drifting  room  for  deputy  satellite  trajectories  and  not  surprisingly  result  in  valid 
trajectories  that  have  a  larger  time-of-flight.  The  overall  shape  of  the  lobe  changes 
little  while  the  peak  rises  rather  substantially  with  larger  and  Ty  (Figure  25).  Note 
that  the  peak  time-of-flight  does  not  increase  linearly  with  lobe  size.  In  our  example 
we  see  a  50%  increase  by  doubling  Tx,  Ty  and  a  200%  increase  by  tripling  them.  Since 
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the  lobe  shapes  do  not  change  appreciably,  we  expect  that  the  optimal  trajectories 
will  be  similar  albeit  with  a  larger  total  time-of-flight.  This  is  shown  in  Section  5.6. 


Figure  25:  Constraint  Surfaces  for  Varying  Size  Lobes  {a  =  45°,  7  =  2  km,  rj  =  45°) 


3.6.7  Sensitivity  of  the  Constraint  Surfaee  to  x^.  Increasing  the  X  coordi¬ 
nate  of  the  lobe  center  {xif)  also  has  an  effect  on  the  constraint  surface.  Recall  that 
a  nonzero  a  value  represents  a  difference  in  the  chief  and  deputy’s  semi-major  axis, 
resulting  in  drifting  orbits.  The  larger  this  difference,  the  greater  the  difference  in 
orbital  period  and  the  bigger  the  relative  velocities.  Thus,  even  though  trajectories 
themselves  appear  similar,  the  relative  velocities  are  larger  in  lobes  further  from  the 
Y  axis,  causing  a  decrease  in  the  maximum  time-of-flight  peak  (Figure  26).  Of  partic¬ 
ular  note,  it  appears  that  doubling  xl  has  the  same  reducing  effect  on  the  constraint 
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Figure  26:  Constraint  Surfaces  for  Varying  Distances  from  the  Y  Axis  = 

1  km,  Ty  =  0.5  km,  rj  =  45°) 


surface  as  halving  the  size  of  the  lobe.  Further  study  is  warranted  and  will  be  left  to 
future  researchers. 


3.6.8  The  Constraint  Surface  for  Elliptical  Chief  Orbits.  Precomputing  the 
constraint  surface  for  elliptical  chief  orbits  is  similar  to  the  circular  case  with  two 
notable  exceptions.  First,  the  elliptical  relative  equations  of  motion  (Equation  (6)) 
are  functions  of  true  anomaly  {p)  and  thus  constraint  surfaces  are  only  valid  for  a 
specific  initial  Po-  This  means  that,  at  the  end  of  a  leg,  a  new  constraint  surface 
with  a  new  Po  must  be  used  to  find  Tmax  for  the  next  leg.  For  a  given  lobe  and 
eccentricity,  the  chief’s  orbit  must  be  discretized  with  appropriate  resolution  and  a 
library  of  constraint  surfaces  stored  in  order  to  hnd  multiple  leg  trajectories.  If  only 
a  single  leg  is  needed,  a  single  constraint  surface  can  be  used. 
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Second,  unlike  the  circular  chief  case,  there  is  no  closed-form  solution  for  the 
relative  velocities.  Several  researchers  (see  Section  2.5  of  the  literature  review)  have 
proposed  high  fidelity  models  for  finding  relative  velocities  for  elliptical  chiefs,  how¬ 
ever,  a  very  simple  and  robust  method  was  used  to  calculate  relative  velocities  in  this 
research.  Using  final  position  error  as  the  cost  function,  Matlab’s®  FSOLVE  algo¬ 
rithm  can  find  the  relative  velocity  given  an  initial/final  position  and  a  time-of-flight 
with  acceptable  accuracy  and  speed. 
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Figure  27:  Elliptical  Chief  Constraint  Surface  e  =  0.7,  Pq  =  Q  rad 
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Figure  28:  Elliptical  Chief  Constraint  Surface  e  =  0.7,  Uq  =  tv  rad 

Notice  the  striking  difference  between  the  constraint  surface  magnitudes  at 
perigee  (Figure  27)  and  apogee  (Figure  28).  This  is  due  to  the  corresponding  dif¬ 
ference  in  relative  velocity  magnitudes  at  apogee  and  perigee.  At  apogee,  the  inertial 
and  relative  velocities  have  a  smaller  magnitude  than  at  perigee  and  therefore  pass 
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between  any  two  points  considerably  more  slowly,  resulting  in  a  larger  Tmax-  This  phe¬ 
nomena  must  factor  into  the  overall  hovering  strategy  since  more  burns  with  higher 
Al/’s  will  be  required  for  hovering  at  perigee  as  opposed  to  apogee. 

3.6.9  Other  Constraints.  Certainly  there  are  other  attractive  constraints 
that  could  be  applied  to  the  cost  function  in  order  to  coax  even  more  desirable  tra¬ 
jectories.  Perhaps  the  first  that  comes  to  mind  is  an  equality  constraint  or  minimum 
constraint  on  total  time-of-flight  in  order  to  satisfy  a  mission  objective,  such  addi¬ 
tions  should  be  avoided.  Any  constraint  added  to  an  optimal  control  problem  risks 
obscuring  a  candidate  extremal  point.  Mission  planners  can  specify  the  number  of 
legs  and  can  control  if  the  trajectory  closes  back  on  itself  (repeating  hover  orbit  con¬ 
dition),  both  of  which  have  the  potential  to  increase  time-of-flight  without  imposing 
additional  constraints.  The  repeating  hover  orbit  in  particular  has  the  potential  to 
satisfy  any  desired  hover  time  by  traversing  the  same  trajectory  over  and  over  again 
until  the  minimum  time-of-flight  is  met. 

3.6.10  Summary.  Examination  of  the  constraint  surface  provides  key  insight 
into  the  behavior  of  the  solution  as  lobe  parameters  change.  Some  are  intuitive  such 
as  larger  lobes  yield  larger  maximum  time-of-flight  and  thus  have  the  potential  to 
lower  the  cost  function.  Others  are  less  intuitive  such  as  how  changes  in  orientation 
angle  which  shift  the  peak  location  will  affect  the  solution.  These  will  be  explored 
in  Chapter  V.  Now  that  we  have  derived  the  equations  of  motion,  developed  a  cost 
function,  and  bounded  the  problem  with  a  constraint,  we  can  start  formulating  a 
trajectory  planner  to  meet  our  stated  objectives.  The  ultimate  goal  is  to  produce 
an  algorithm  that  outputs  a  set  of  impulsive  thrust  locations  and  the  time-of-flight 
between  them.  This  is  done  in  the  next  chapter. 
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IV.  The  Optimal  Trajectory 

The  optimal  trajectory  will  be  the  output  of  a  nonlinear  programming  algo¬ 
rithm  that  searches  for  the  lowest  value  of  the  cost  function  while  staying  within 
the  maximum  time-of-flight  constraint.  We  seek  a  method  in  which  the  user  pro¬ 
vides  the  number  of  legs  (k),  entry/exit  conditions,  and  lobe  position,  orientation, 
and  shape,  after  which  the  algorithm  provides  the  angular  positions  of  the  impulsive 
thrust  locations  and  the  time-of-flight  between  them: 

Given  Output 

k  (Number  of  legs)  ,  •02 ,  •  •  • ,  t/'fc+i 

Entry  Condition  f'1,2,  T2,3,  ...,fk,k+i 

Exit  Condition 
a,  p,  7 
Tc)  T/’  ^ 

These  outputs  can  easily  be  transformed  into  impulsive  thrust  vectors  and  thrust 
times.  The  propagated  trajectories  shown  in  the  results  chapter  are  calculated  by 
taking  the  hnal  velocity  of  the  previous  leg  (or  the  entry  conditions  for  the  first  leg), 
adding  the  impulsive  thrust  vector  and  then  propagating  until  the  next  thrust  time. 
In  this  manner,  we  can  ascertain  the  error  incurred  by  using  a  reduced  order  model 
for  the  optimization  and  then  applying  it  to  a  higher  fidelity  truth  model. 

Since  motion  in  the  XY  plane  and  Z  direction  decouple,  the  optimal  trajectories 
for  both  will  be  found  separately.  Optimization  in  the  XY  plane  produces  a  trajectory 
and  total  time-of-flight  (Tr).  This  total  time-of-flight  is  then  used  to  find  a  closed- 
form  optimal  trajectory  for  motion  in  the  Z  direction.  The  two  motions  can  then  be 
combined  for  a  final  trajectory.  We  start  with  optimization  in  the  XY  plane. 
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4-1  The  Optimal  Trajectory  in  the  XY  Plane 

The  cost  function  is  Equation  (20)  from  Section  3.3: 


oa,j.2  jr-2p2 

j  _  OO-^min'*  -‘-T 


AE2  +  ^y2  ^  ^y2  ^  ^  ^y2  ^  ^y2 


Tio  +  ^2  3  +  T3  4  + ...  + 


Tp 


where  AV  is  specific  AV  and  has  no  dependence  on  the  chief’s  semi-major  axis  or  fi. 
In  an  effort  to  keep  the  numerator  close  to  unity,  the  AE’s  are  scaled  by  the  amount 
of  specific  AV  required  to  remain  at  the  minimum  x  coordinate  of  the  lobe  boundary 
(3:^min)-  If  the  lobe  intersects  the  Y  axis  (xmin  =  0),  an  alternative  optimal  trajectory 
is  employed  and  is  discussed  in  Section  5.1.  The  constraint  on  each  time-of-fiight  is 
(from  Section  3.6) 

Ti,i+1  A  (V’l ;  ^^2) 

where  T^ax  is  the  maximum  time-of-fiight  that  keeps  the  entire  trajectory  within  the 
lobe.  There  are  no  constraints  on  the  angular  positions  ('^)  unless  specified  by  the 
entry/exit  conditions  which  are  discussed  later  in  this  chapter.  The  cost  function  is 
minimized  via  the  FMINCON  nonlinear  programming  routine  in  Matlab®  which 
is  specifically  designed  for  nonlinear  cost  functions  that  have  nonlinear  constraints. 
The  pseudocode  for  this  routine  is  found  in  Appendix  Q. 


4-2  The  Optimal  Trajectory  in  the  Z  Direction 

Recall  that  motion  in  the  Z  direction  is  a  harmonic  oscillator  with  a  closed-form 
solution  of  (Equation  (9),  Section  3.2) 

'^(t)  '^max  cos  (27rt  ~\~  0) 


where 


=  tan 


-1 


2t^Zo 
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It  can  be  shown  (derivation  in  Appendix  J)  that  the  maximum  time-of-flight  in  which 
the  deputy  stays  between  a  minimum  (5min)  and  maximum  (5max)  Z  coordinate  is 

f„„  =  ieos-‘(|=-')  (26) 

^  \  ^max  / 

The  cos“^  function  has  a  domain  of  0  <  I™  <  1  and  this  fits  well  with  our  definition 
of  the  elliptical  cylinder.  Motion  in  the  Z  direction  is  caused  by  differences  in  the 
chief  and  deputy’s  inclinations  and/or  longitude  of  the  ascending  node.  As  intuition 
and  Equation  (25)  show  us,  the  period  of  the  Z  oscillation  is  27r  and  is  centered  about 
the  chief’s  orbit  plane.  Half  of  the  chief’s  period  is  spent  above  the  orbit  plane  and 
half  below.  Thus,  any  Z  period  greater  than  0.5  would  indicate  that  the  deputy  has 
passed  through  the  chief’s  orbit  plane.  Lobes  that  include  the  chief’s  orbit  plane  will 
have  .?min/5max  ratios  that  are  less  than  zero  (change  in  sign  between  z^iin  and 
and  ill-defined  lobes  are  greater  than  1  (|zrnin|  is  larger  than  |5max|)-  Therefore,  either 
the  lobe  intersects  the  orbit  plane,  in  which  case  the  optimal  solution  is  to  stay  in  the 
orbit  plane  with  no  Z  motion,  or  the  lobe  is  improperly  defined.  We  can  expect  only 
ratios  between  0  and  1  and  Tmax’s  between  0  and  0.5  chief  orbits. 

A  harmonic  oscillator  affords  a  single  degree-of-freedom  to  optimize  the  problem. 
The  period  of  the  Z  oscillation  (P^)  is  chosen  as  the  optimization  parameter  based 
on  the  following  argument.  Assume  the  XY  plane  optimization  has  yielded  a  total 
time-of-flight  (2t)-  We  desire  the  deputy  to  stay  between  z^iin  and  z^ax  for  the  same 
amount  of  time.  The  minimum  number  of  burns  required  is 


^  Burns 


Pz 


where  [  J  represents  the  floor  function.  It  can  be  shown  that  the  AV  required  for 
each  of  these  burns  is  (derivation  in  Appendix  K.3) 


AI/  = 


AE 


n 


=  25min  tan  (  ttP 
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Thus  the  total  AV  required  is 


AVz 


,  tan  I  ttP, 


25min  {kz 


(27) 


where  kz  is  the  number  of  legs  in  the  Z  direction.  Note  that  AVz  is  inversely  propor¬ 
tional  to  Pz\  thus  we  desire  the  largest  possible  period  (or  fewest  number  of  bounces, 
kz)  without  exceeding  the  maximum  period  dehned  by  Equation  (26).  An  example 
problem  is  found  in  Figure  29.  Assume  that  the  XY  optimization  has  yielded  a  total 
time-of-flight  of  Tt  =  0.45,  that  =  1,  and  =  1.25.  The  period  of  the  Z 
motion  is  therefore  constrained  to  be  less  than 


—cos  ^  (0.8)  ~  0.2 

TT 


Referencing  Figure  29,  the  smallest  possible  total  AV^  occurs  at 


Tt 

P,  =  —  =  0.15 
3 


The  green  line  is  the  AV  required  for  a  continuous  thrust  hover  at  and  is  equal 
to  (Equation  (137)  in  Appendix  L) 


AVc  = 


Note  how  the  discrete  solution  converges  to  the  continuous-thrust  line  as  0, 

that  is,  as  the  number  of  bounces  (kz)  goes  to  infinity  (see  Appendix  K.3  for  the  full 
proof). 
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Figure  29:  Notional  AV  vs  Z  Period 

4-3  Entry  Condition  Definitions 

In  addition  to  lobe  size  and  orientation,  the  mission  planner  has  flexibility  in 
choosing  entry  and  exit  conditions  to  and  from  the  lobe.  Combinations  of  these 
conditions  are  used  to  define  the  cases  presented  later.  There  are  three  types  of  lobe 
entry  conditions  used  in  this  research.  They  are: 

1.  Defined  Entry  Condition 

2.  Open  Entry  Condition 

3.  Entry  from  a  Closed  Relative  Orbit 
and  are  explained  in  the  following  subsections. 

4.3.1  Defined  Entry  Condition.  The  first  entry  condition,  Defined  Entry 
Condition  (DEnC),  enables  the  user  to  specify  the  deputy’s  relative  entry  position 
and  velocity.  Thus,  in  addition  to  the  standard  inputs  listed  at  the  beginning  of  this 
chapter,  the  user  must  specify: 
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This  will  affect  the  calculation  of  AV^  (Equation  (124)  in  Appendix  K): 


_  ry  n 

1  Ul 


Vi 


00  00 

yt  -  Vi 


where 


A  + 

— 

Vi 

-yf_ 

-45++67rTi,2C+ 

8-67rTi,25+-8C+ 

-14+127rfi,2.g++14C+ 

8-67rTi,25+-8C+ 


45+-67rf'i,2 

-67rTi,25+-8C+ 

2-2(7+ 


-2+2(7+ 


-67rTi,2S'+-8C+ 

S+ 


8-67rTi,2S'+-8C+  8-67rTi,2S'+-8C+ 


1 

Xi 

r  1 

Xi 

X2 

— 

Vi 

Ay 

Ay  =  ^2  -  yi 


4. .3. 2  Open  Entry  Condition.  The  second  variant,  Open  Entry  Condition 
(OEnC),  pnts  no  restrictions  on  the  relative  start  position  and  velocity;  thus  the 
algorithm  is  free  to  choose  these  in  whatever  manner  is  most  optimal.  Since  the 
deputy  is  entering  the  lobe  at  the  desired  position  and  velocity,  there  is  no  thrust  at 
the  first  point: 

=  0 


4.3.3  Entry  from  a  Closed- Relative  Orbit.  The  Entry  from  a  Closed- Relative 
Orbit  (EnCRO),  enables  the  algorithm  to  search  for  an  optimal  start  position  but 
constrains  the  velocity  at  that  point  to  match  a  closed-relative  orbit  that  is  centered 
on  the  chief  and  intersects  the  lobe  at  that  point.  While  closed-relative  orbits  can 
be  centered  about  any  point  along  the  Y  axis,  the  choice  was  made  to  center  on  the 
chief.  Setting  the  relative  orbit  parameters  a  and  b  (Equations  (10a)  and  (10b)  in 
Section  3.2)  to  zero. 


xi  =  Tiyi 
y-  =  -YkXi 


56 


This  will  effect  the  calculation  of  AV^: 


_  ry  rt  i~^ 

1  iJl 


00  00  ^ 

yt  -  vO 


where 


xf 

—  Xi 

-4S'++67rTi,2C'+ 

45+-67rri,2 

-2+2(7+ 

8-67rTi,25+-8C+ 

8-67rTi,25+-8(7+ 

8-67rTi,2  5+-8(7+ 

-yi  _ 

- 14+ 1 27rf'i  ,2  S+ + 14(7+ 

2-2(7+ 

5+ 

.  8-67rTi,25+-8(7+ 

8-67rTi,25+-8(7+ 

8-67rTi,2  5+-8(7+. 

TT^l 

—Anxi 


^y  =  y2-  m 


4-4  Exit  Condition  Definitions 

The  lobe  exit  conditions  nsed  in  this  research  are 

1.  Defined  Exit  Condition 

2.  Open  Exit  Condition 

3.  Exit  to  a  Closed  Relative  Orbit 

4.  Repeating  Hover  Orbit 

and  are  explained  in  the  following  subsections. 

4-4-i  Defined  Exit  Condition.  The  first  exit  condition,  Defined  Exit  Con¬ 
dition  (DExC),  enables  the  nser  to  specify  the  depnty’s  relative  exit  position  and 
velocity.  Thns,  in  addition  to  the  standard  inpnts  listed  at  the  beginning  of  this 
chapter,  the  user  must  specify: 


^Pk+i, 
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This  will  effect  the  calculation  of  ^Vp  (Equation  (124)  in  Appendix  K): 


= 


x{  -x^  Vi  -y^ 


-1 

Ry  1  Ry  ^ 

•  1 

y{  -  Vi  _ 

where 


xf 

—  x-^ 

rp  1 

yt 

-yi  _ 

yt  _ 
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45+-67rfi,2 

-2+2(7+ 

Xi 

8-67rfi,2S+-8C'+ 

-14+127rTi,25++14C+ 

8-67rfi,2S+-8C'+ 

2-2C+ 

8-67rfi,2S+-8(7+ 

S+ 

X2 

8-67rfi,25+-8C+ 

8-67rTi,25+-8C+ 

8-67rfi,25+-8(7+. 

Ay 

Ay  =  ^2  -  m 


4.4-2  Open  Exit  Condition.  The  second  variant,  Open  Exit  Condition 
(OExC),  puts  no  restrictions  on  the  relative  exit  position  and  velocity.  Thus  the 
deputy  will  leave  the  lobe  on  the  trajectory  resulting  from  the  AVk  burn.  There  is  no 
burn  at  the  final  point: 

AV^  =  0 

4-4-3  Exit  to  a  Closed-Relative  Orbit.  The  Exit  to  a  Closed- Relative  Orbit 
(ExCRO),  enables  the  algorithm  to  search  for  an  optimal  exit  position  but  constrains 
the  velocity  at  that  point  to  match  a  closed-relative  orbit  that  is  centered  on  the  chief 
and  intersects  the  lobe  at  that  point.  While  closed-relative  orbits  can  be  centered 
about  any  point  along  the  Y  axis,  the  choice  was  made  to  center  on  the  chief.  Setting 
the  relative  orbit  parameters  a  and  b  (Equations  (10a)  and  (10b)  in  Section  3.2)  to 
zero. 


^k+l  T^Vk+l 

Vk+I  =  -4™^+! 
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AV^  is: 


AV^  = 


^k+l  ^k+1  Vk+l  Vk+l 


^k+1  ^k+1 

^t+1  -  ^k+1 


where 


'^k+l 

^k+l 

-45-+67rf 

45--67rfC- 

2-2(7- 

_ 

8-67rf5--8C- 

8-67rT5--8C'- 

8-6nfS--8C- 

1 - 

+  + 

-  ^k+l_ 

-^TTXk+l 

to 

1 

to 

1 

-14+127rfS-+14(7 

s- 

.8-6-KfS--8C- 

8-67rT5--8C- 

8-6nfS--8C- _ 

Xk 

Xk+l 

Ay 


Ay  =  Vk+l  -  Vk 


4-4-4  Repeating  Hover  Orbit.  The  repeating  hover  orbit  is  not  an  exit 
condition  per  se.  Instead,  we  seek  a  closed  relative  orbit  contained  within  the  lobe  in 
which  the  starting  relative  position  and  velocity  are  identical  to  the  ending  conditions. 
This  prodnces  a  trajectory  that  can  be  traversed  for  as  many  periods  as  desired  or  as 
long  as  fuel  is  available.  In  order  to  enforce  a  closed  orbit,  the  following  constraints 
are  applied  to  the  cost  function: 


Xk+l  X\ 
Vk+l  =  yi 
Tk+i,k  =  0 
Tp  =  Tk,i 


Xk 

Vk 


AV^  = 


Xk  Vk  xi  yi  X2  2/2 


Xi 

Vi 


X2 


y2 
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4^5  Research  Cases 

The  availability  of  three  entry  and  four  exit  conditions  provides  us  with  twelve 
possible  entry /exit  combinations.  In  the  interest  of  brevity  and  given  that  only  a 
handful  of  combinations  are  of  true  interest  to  mission  planners,  only  three  cases  are 
presented  herein.  They  are: 

1.  Defined  Entry  Condition/Open  Exit  Condition 

2.  Entry  From  a  Closed  Relative  Orbit/Exit  to  a  Closed  Relative  Orbit 

3.  Open  Entry  Condition/Repeating  Hover  Orbit 

4' 6  The  Continuous- Thrust  Solution 

As  a  benchmark  to  which  to  compare  our  discrete-thrust  trajectories,  a  continuous- 
thrust  solution  is  derived.  We  can  quickly  find  the  closed-form  solution  for  a  continuous- 
thrust  controller  that  keeps  the  deputy  at  a  specified  point  in 
shortened  version  of  the  derivation  in  Appendix  L  is  provided 
at  this  hover  point  must  be  zero: 

X  =  0  y  =  0  z  =  0 

which  reduces  Equation  (7)  to: 

X  =  47r(0)  -|-  127r^5:  =  127r^7: 
y  =  — 47r(0)  =  0 
5  = 

Since  the  y  equation  is  now  zero,  we  need  only  worry  about  the  accelerations  along 
the  X  and  Z  axes.  Integrating  from  zero  to  the  total  time-of-flight  (TV)  yields  the 
AV  as  a  function  of  chief  orbit  period  (Al^)  required  to  keep  the  deputy  hovering  at 


the  relative  frame.  A 
below.  The  velocities 
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a  given  yo,  Zo)\ 


.Tt 


r^T 


cTt 


rTr 


AVc  = 


\x\dt  +  /  \z\dt  =  /  12n^\xo\dt+  /  4:7T^\zo\dt 


rTr 


=  {l2\Xo\+  ^\Zo\)  /  dt  =  {l2\Xo\+ ^\Zo\)Tt 


where  the  C  subscript  indicates  “continuous” .  The  location  of  Xo  and  Zo  is  arbitrary 
but  the  smallest  continuous  Al/  is  attained  when  Xo  and  Zo  are  at  their  minimum, 
thus  5;min  and  z^iin  represent  the  absolute  value  of  the  coordinate  on  the  lobe  that  is 
closest  to  the  Y  axis  and  XY  plane  respectively,  and 


Al/c  =  +  45i„in)  TT^TV 


(28) 


If  the  lobe  happens  to  intersect  the  Y  axis,  then  the  optimal  solution  (for  the  linear 
CW  case)  is  to  stay  on  that  Y  axis  (the  loci  of  eqnilibrinm  points)  or  on  a  closed  2x1 
ellipse  about  it  (see  Section  5.1).  We  can  also  express  the  continnous  Al^  as  a  specific 
AV  by  applying  the  conversion  in  Equation  (17),  Appendix  3.3 

(12Xmin  +  4Zmin)  =  {Qx^iu  +  S^min)  TtTt  (29) 

The  above  eqnation  works  well  if  we  assnme  the  deputy  starts  at  the  minimum  x  and 
values,  however,  in  order  to  use  this  as  a  fair  comparison,  the  AV  needed  to  get  into 
and  ont  of  that  position  mnst  also  be  inclnded.  This  motivates  finding  minimum- fuel 
entry  and  exit  legs  as  shown  in  Fignre  30. 


Al/c  = 


AVc 

27r 


4.6.1  Entry  Leg.  Assume  that  the  continuous  burn  solution  starts  at  the 
same  entry  position  (xi,yi)  and  velocity  {xi,yi)  as  the  discrete-thrust  solution.  The 
AV  required  to  get  to  the  Xmin  position  is  (from  Equation  (124)  in  Appendix  K): 
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Figure  30:  The  Continuous-Thrust  Trajectory 
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Once  at  the  x^in  position,  the  deputy  must  make  a  burn  to  cancel  its  relative  velocity 
that  is  equal  in  magnitude  but  opposite  in  direction  to  the  final  velocity  of  the  previous 
leg 

AV2  =  ^\\i2X  +  ^2n2 
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We  can  use  either  a  minimization  routine  or  a  global  search  over  the  range  of  time- 
of-flight: 

0  <  Ti,2  <  7maxi,2 

to  find  the  smallest  total  AV  required  to  arrive  and  stop  at  x^m'- 


A^Entry  =  min 


AW  +  AW 


4.6.2  Exit  Leg.  The  exit  leg  will  depend  on  the  user-specified  exit  condition. 
If  the  exit  condition  is  open,  the  deputy  is  allowed  to  drift  out  of  the  lobe  and  requires 
no  additional  AW : 

AWxit  =  0 

If  a  repeating  hover  orbit  is  requested,  the  deputy  never  leaves  the  Xmin  position  and 
AWxit  is  again  zero.  If  the  exit  condition  is  an  exit  to  a  closed-relative  orbit  or  if  the 
exit  position  and  relative  velocity  are  specified,  then  the  deputy  will  need  to  perform 
two  burns  to  leave  x^m  to  arrive  at  the  exit  position  and  then  to  leave  the  lobe  at  ■03 . 
The  AW  required  to  get  from  the  x^in  position  to  the  exit  point  is  (from  Equation 
(124)  in  Appendix  K) 

Aq  =  i|iq.t  +  ^+y||2 
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The  AW  required  to  exit  the  lobe, 
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where  for  the  defined  exit  condition 
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For  the  exit  to  a  closed-relative  orbit,  the  relative  velocity  at  ips  must  be  the  same  as 
a  closed- relative  orbit  centered  on  the  chief  that  passes  through  Setting  a  =  b  =  0 
(Equations  (10a)  and  (10b)), 
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We  can  use  either  a  minimization  routine  or  a  global  search  over  the  range  of  time- 
of-flight: 

0  <  T2^3  <  T’max2,3 


to  find  the  smallest  total  AV  required  to  leave  Xmin  and  exit  the  lobe: 


ACexh  =  niin 


ACs  +  AW 


4-6.3  Hover  Time.  The  time-of-fiight  for  each  leg  is  obviously  independent 
of  the  total  time-of-flight  produced  by  the  discrete-thrust  algorithm.  Whatever  time 
remains  after  getting  to  Xmin  and  exiting  the  lobe  will  be  spent  hovering  at  Xmin  with 
continuous-thrust;  thus: 


^Hover 


—  Tt  —  ^1,2  —  ^2^3 
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There  are  cases  in  which  the  discrete-thrust  solution  total  time-of-flight  is  less  than 
the  sum  of  the  two  continuous-thrust  legs,  this  must  be  taken  account  when  comparing 
the  two  results. 

4.6.4  Summary.  This  chapter  has  derived  all  the  necessary  components  to 
generate  fuel-optimal  discrete-thrust  trajectories  as  well  as  a  benchmark  continuous 
method  with  which  to  compare  those  results.  We  are  finally  ready  to  present  and 
analyze  those  final  results. 
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V.  Results  and  Analysis 


With  an  understanding  of  the  dynamics  and  solution  method  in  hand,  we  turn 
next  to  analyzing  the  results  and  drawing  conclusions.  The  symmetries  of  the  con¬ 
straint  and  cost  function  assure  us  that  optimal  trajectories  for  lobes  in  the  upper 
right  quadrant  are  representative  of  lobes  in  any  other  quadrant  (see  Section  3.6). 
This  is  very  convenient,  as  it  allows  us  to  use  a  small  number  of  results  to  make 
general  conclusions  about  optimal  trajectories.  Further,  these  results  are  applicable 
to  any  sized  circular  chief  orbit  and  about  any  celestial  body,  assuming  two  body 
dynamics  are  the  dominant  force. 

First,  results  are  presented  for  the  a  special  class  of  hovering  orbits  that  take 
advantage  of  the  equilibrium  condition  of  the  equations  of  motion.  Next  are  optimal 
trajectories  for  each  of  the  three  cases  defined  in  Section  4.5  as  applied  to  lobes  2, 
3,  and  4  (shown  in  Figure  31  with  parameters  in  Table  1).  Lobe  1  is  a  special  case 
in  which  the  lobe  intersects  the  Y  axis  and  will  only  be  examined  in  the  persistent 
hover  orbit  section  (Section  5.1).  Trajectories  are  propagated  with  the  unperturbed, 
linearized  equations  of  motion. 


''' . • 

. -J 

'  ''''  / 

/  /  \ 

1 

1 

1 

1 

1 

1 

1 

f 

-1000  -500  0  500  1000  1500  2000  2500  3000 

X  Axis  (m) 


Figure  31:  Lobes  of  the  Results  Section 
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Table  1:  Lobe  Parameters 


Lobe 

# 

a 

(degs) 

(degs) 

7 

(km) 

Tx 

(km) 

(km) 

V 

(degs) 

1 

90 

90 

2 

1 

0.5 

20 

2 

45 

90 

2 

1 

0.5 

45 

3 

45 

90 

2 

0.5 

1 

45 

4 

0 

90 

1.414 

1 

0.5 

0 

5.1  The  Persistent  Hover  Orbit 

As  mentioned  earlier  and  discussed  in  detail  in  Section  3.2.2  and  Appendix  M, 
the  Y  axis  is  the  locus  of  equilibrium  points  for  the  linearized  equations  of  motion 
about  a  chief  in  a  circular  orbit  (Equation  (7)).  Lobes  that  contain  some  portion 
of  the  Y  axis  deserve  special  consideration.  First,  let  us  examine  the  output  of  the 
optimization  algorithm  when  given  a  lobe  that  intersects  the  Y  axis.  Figure  32  is  the 
result  of  an  optimization  run  for  a  defined  entry  {'ipi  =  |)  and  open  exit  on  Lobe  1. 
The  entry  relative  velocity  is  chosen  such  that  the  deputy  enters  from  a  closed  relative 
orbit  centered  on  the  chief  at  ipi.  The  result  is,  not  surprisingly,  to  place  the  deputy  on 


Figure  32:  Optimization  Algorithm  Results  for  Lobe  1 

a  2x1  ellipse.  This  special  case  relative  orbit  is  called  a  persistent  hover  orbit  (PHO), 
defined  as  any  natural  closed  orbit  that  fits  entirely  within  the  lobe  (see  Figure  33). 
The  PHO  is  a  type  of  repeating  hover  orbit  that  does  not  require  any  impulsive 
thrusting  to  maintain,  other  than  what  is  required  to  correct  for  perturbations  and 
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Figure  33:  The  Persistent  Hover  Orbit 

linearization  errors,  making  it  ideal  for  long-term  proximity  operations.  The  PHO  is 
specified  by  two  parameters:  the  semi-minor  axis  (p)  and  the  center  coordinate  along 
the  Y  axis  (6).  In  the  limit  as  p  — 0,  the  PHO  becomes  a  point  on  the  Y  axis.  Finding 
acceptable  PHOs  for  a  given  lobe  is  the  subject  of  Appendix  M.  The  PHO  problem  is 
as  follows.  Assume  the  deputy’s  trajectory  intersects  the  lobe  at  some  entry  position 
with  a  given  relative  velocity.  Positioning  the  deputy  onto  the  PHO  requires  two 
impulsive  thrusts,  the  first  at  lobe  entry  to  place  the  deputy  on  a  trajectory  towards 
the  PHO  and  another  to  enter  it.  This  greatly  simplifies  our  optimization  routine 
since  we  can  globally  search  for  the  minimum- fuel  solution  of  these  two  burns.  The 
only  remaining  question,  is  what  are  the  best  values  of  the  PHO  semi-minor  axis  (p) 
and  center  position  (6). 

Although  an  exhaustive  proof  will  not  be  made  here,  indications  are  that  the 
minimum- fuel  answer  is  to  choose  p  and  b  such  that  p  is  as  large  as  possible.  Figures 
34  and  35  present  the  result  of  evaluating  all  values  of  p  and  b  that  yield  PHOs  that 
are  completely  contained  inside  Lobe  1.  The  left-hand  side  figures  are  the  AV 
surface  over  p  and  b  while  the  right-hand  figure  shows  the  resultant  minimum-fuel 
trajectory.  Note  that  only  combinations  of  p  and  b  that  fit  completely  inside  the 
lobe  are  evaluated  and  the  stairstep  feature  is  due  to  discretization  of  the  search 
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(b)  The  Minimum  AV  Trajectory 


Figure  34:  Persistent  Hover  Orbit  Test:  j 


(a)  AF  as  a  Function  of  p  and  b 


(b)  The  Minimum  AV  Trajectory 


Figure  35:  Persistent  Hover  Orbit  Test:  ipi  =  ^ 


space.  In  both  cases,  the  minimum  AV  occurred  when  the  largest  possible  p  value 
was  chosen.  Note  that  if  mission  requirements  dictate  a  fairly  small  time-of-flight, 
then  optimizing  the  cost  function  over  one  or  two  legs  may  still  yield  the  lowest  total 
fuel  solution;  however,  if  a  longer  total  time-of-flight  is  required,  searching  for  the 
minimum  two-thrust  solution  to  place  the  deputy  onto  a  PHO  is  optimal. 
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5.2  Case  1:  Defined  Entry  Condition/ Open  Exit  Condition 

Case  1  investigates  trajectories  for  which  the  entry  condition  is  defined  and  the 
exit  condition  is  left  open.  Recall  that  the  defined  entry  condition  requires  the  user  to 
provide  an  incoming  relative  velocity  at  a  specified  lobe  entry  point.  For  convenience, 
and  because  choice  of  this  initial  condition  is  arbitrary,  the  following  results  will  use 
the  relative  velocity  of  the  closed  relative  orbit  that  intersects  the  lobe  at  the  specified 
point  and  is  assumed  to  be  centered  on  the  chief.  Setting  the  relative  orbit  parameters 
a  and  b  (Equations  (10a)  and  (10b),  Section  3.2)  to  zero  yields 

y-  =  -AtTXi  Xjj  =  TT^i 

Case  1  is  of  particular  interest  because,  for  the  single-leg  case,  there  are  only  two 
optimization  parameters  (^2  and  Tj,2);  thus  we  can  plot  the  cost  surface  and  watch 
the  iteration  path  of  the  nonlinear  programming  algorithm.  Let 


The  results  of  the  optimization  algorithm  for  each  of  the  three  lobes  under  the  con¬ 
ditions  specified  above  are  shown  in  Figures  36,  37,  and  38  where  the  red  line  on  the 
cost  contour  plot  is  the  line  of  maximum  time-of- flight. 

In  all  three  examples,  we  can  visually  confirm  that  the  optimization  routine  has 
converged  on  a  minimum  solution.  We  can  also  confirm  that  these  are  not  in  fact 
the  global  minima  for  this  particular  entry/exit  condition.  In  each  lobe,  the  closed 
relative  orbit  on  which  the  deputy  arrives  can  pass  through  the  lobe  and  intersect  a 
second  point  on  the  lobe.  Since  the  exit  condition  is  open,  the  zero-fuel  solution  is  to 
stay  on  the  original  closed  relative  orbit  and  exit  the  lobe  at  the  second  intersection 
point.  This  is  most  clearly  demonstrated  in  Lobe  3  (Figure  37)  for  which  the  true 
global  minimum  is  'ip2  =  5.31  rad,  Ti^2  =  0.048.  Figure  39  is  a  repeated  run  of  the 
optimization  algorithm  with  different  initial  conditions  that  converges  to  the  global 
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(a)  Cost  Function  and  Optimization  Itera¬ 
tions 


Figure  36:  Case  1,  Lobe  2,  Single  Leg 
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(a)  Cost  Function  and  Optimization  Itera¬ 
tions 


(b)  Optimal  Discrete  Trajectory 


Figure  37:  Case  1,  Lobe  3,  Single  Leg 

minimum  above.  This  is  the  trivial  solution  and  is  not,  in  general,  the  solution  of 
interest  for  the  mission  planner. 

The  example  above  highlights  the  sensitivity  of  the  optimization  algorithm  to 
initial  condition  choice.  As  with  most  nonlinear  optimization  problems,  we  must 
be  cognizant  of  the  possibility  of  multiple  local  minima.  Two  are  apparent  in  our 
simple  example  above  and  more  should  be  expected  as  the  dimensionality  increases 
(by  increasing  the  number  of  legs,  k).  Unfortunately,  we  do  not  have  the  luxury  of 
visually  verifying  convergence  to  the  correct  solution  when  going  to  higher  dimensional 
search  spaces.  In  order  to  increase  confidence  in  finding  the  desired  minimum  point. 
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(a)  Cost  Function  and  Optimization  Itera¬ 
tions 


Figure  38:  Case  1,  Lobe  4,  Single  Leg 
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(a)  Cost  Function  and  Optimization  Itera¬ 
tions 


(b)  Optimal  Discrete  Trajectory 


Figure  39:  Case  1,  Lobe  3,  The  Trivial  Solution 


all  simulations  were  execnted  with  mnltiple  initial  guesses  with  the  lowest  cost  result 
presented.  Currently  the  initial  choice  of  xp  is  accomplished  by  the  user;  however, 
there  is  no  reason  the  process  could  not  be  automated  in  order  to  look  over  a  wide 
range  of  initial  conditions. 

The  next  set  of  graphs  show  the  evolution  of  trajectories  for  each  of  the  three 
lobes  as  the  number  of  legs  increases.  The  left-hand  plots  show  the  lowest  cost  discrete 
trajectory  found  after  multiple  initial  condition  choices.  The  right-hand  plots  are  the 
time  history  of  the  discrete-  and  continnous-thrust  trajectories.  The  continuous- 
thrust  solution  is  described  in  Section  4.6  and  is  a  function  only  of  the  lobe  shape, 
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position,  orientation,  and  entry  position  and  velocity.  In  general,  it  consists  of  an 
entry  leg  to  the  Xmin  coordinate,  possible  hovering  at  Xmin,  and  exiting  the  lobe  either 
at  a  designated  point,  drifting  out  of  the  lobe,  exiting  onto  another  closed  relative 
orbit,  or  remaining  at  Xmin-  These  thrusts  are  based  on  minimum  AV  solutions.  The 
continuous  solution  entry  and  exit  legs  will  both  have  a  fixed  time-of-flight,  again 
based  on  the  minimum  AV  solution.  The  difference  between  the  sum  of  these  two  legs 
and  the  total  time-of-fiight  of  the  discrete  thrust  solution  is  made  up  by  hovering  with 
continuous-thrust  at  x^m-  Since  there  is  no  minimum  total  time-of-fiight  constraint 
applied  to  the  cost  function,  it  is  possible  that  the  sum  of  the  entry  and  exit  legs 
is  greater  than  the  time-of-fiight  found  in  the  optimization  routine.  In  this  case,  no 
hovering  is  required  and  the  continuous  thrust  time-of-fiight  will  be  longer  than  the 
discrete  solution.  Care  must  be  taken  when  comparing  the  discrete  and  continuous 
solutions  in  these  cases. 

Figures  40-44  show  the  evolution  of  trajectories  for  Lobe  2  as  the  number  of  legs 
increases  from  one  to  five.  There  are  minor  changes  in  the  trajectory  and  performance 
as  additional  legs  are  requested,  but  there  is  a  definite  pattern  of  convergence  towards 
a  final  trajectory.  In  all  but  the  single- leg  case,  the  continuous  solution  outperforms 
the  discrete  case.  Also,  all  but  the  five-leg  simulation  runs  into  the  problem  discussed 
above  with  the  time-of-fiight  of  the  entry  and  exit  legs  of  the  continuous  solution 
being  larger  than  the  entire  discrete  time-of-fiight.  It  is  clear  that,  for  Lobe  1,  a 
continuous-thrust  strategy  is  more  efficient. 

Figures  45-49  are  trajectories  within  Lobe  3.  Again  there  is  clear  convergence 
to  a  final  trajectory,  but  unlike  in  the  previous  case,  there  are  slight  fuel  savings  by 
using  discrete  solution  in  all  but  the  final  [k  =  5)  simulation.  This  set  of  runs  is  also  a 
good  example  of  how  additional  legs  are  accommodated  by  clustering  several  smaller 
thrusts  to  take  the  place  of  one  larger  thrust.  In  most  cases,  this  is  significantly  less 
efficient,  a  fact  notably  displayed  in  Figure  49.  In  the  A:  <  4  trajectories,  the  first  leg 
from  the  entry  point  to  near  the  top  of  the  lobe  is  the  result  of  a  single  thrust.  When 
k  =  5,  this  single  thrust  is  split  into  three  thrusts.  Although  the  time-of-fiight  and 
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Figure  40:  Case  1,  Lobe  2,  Legs  =  1 
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Figure  41:  Case  1,  Lobe  2,  Legs  =  2 


trajectory  of  this  leg  remain  nearly  identical,  the  fuel  spent  is  nearly  70%  greater  in 
the  latter  case.  This  splitting  phenomena  is  due  to  the  unrestricted  bounds  on  each 
leg’s  time-of-flight  and  AC.  With  no  minimum  time-of-flight,  optimal  trajectories 
are  split  into  smaller  pieces  in  order  to  satisfy  the  requirement  for  more  legs.  This 
is  a  shortcoming  of  the  cost  function  as  posed  and  should  be  addressed  in  future 
work.  The  splitting  effect  makes  it  clear  that,  although  there  are  some  increases  to 
time-of-flight,  there  is  a  distinct  k  after  which  there  is  no  improvement. 

The  final  set  of  plots  is  for  Lobe  4.  The  results  are  similar  to  those  above, 
with  the  single  leg  requiring  no  hover  from  the  continuous  solution  and  some  gains  in 
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Figure  42:  Case  1,  Lobe  2,  Legs  =  3 
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Figure  43:  Case  1,  Lobe  2,  Legs  =  4 


time-of-flight  for  higher  numbers  of  legs.  This  scenario  seems  to  split  the  difference 
between  the  hrst  two  as  far  as  performance  comparisons  between  the  continuous  and 
discrete  solutions.  Whereas  the  continuous  solution  provided  lower  fuel  costs  in  Lobe 
2  and  the  discrete  solution  marginally  lower  fuel  costs  in  Lobe  3,  there  are  examples 
of  both  in  Lobe  4.  For  k  <  3,  the  discrete  thrust  solution  requires  slightly  less  fuel 
while  anything  larger  should  use  a  continuous  solution.  The  results  of  Case  1  over  the 
three  lobe  shapes  indicate  that  the  orientation  angle  (rj)  strongly  influences  correct 
method  choice  (discrete  or  continuous). 
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Case  1,  Lobe  2,  Legs  =  5 
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Figure  45: 


Case  1,  Lobe  3,  Legs  =  1 
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Figure  46: 


Case  1,  Lobe  3,  Legs  =  2 
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Figure  47: 


Case  1,  Lobe  3,  Legs  =  3 
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Figure  48: 


Case  1,  Lobe  3,  Legs  =  4 
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Figure  49:  Case  1,  Lobe  3,  Legs  =  5 
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Figure  50:  Case  1,  Lobe  4,  Legs  =  1 
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Figure  52:  Case  1,  Lobe  4,  Legs  =  3 
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Figure  53:  Case  1,  Lobe  4,  Legs  =  4 
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Figure  54:  Case  1,  Lobe  4,  Legs  =  5 
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5.3  Case  2:  Entry  From  a  Closed- Relative  Orbit/ Exit  to  a  Closed- 
Relative  Orbit 

The  entry  from  and  exit  to  a  closed-relative  orbit  add  additional  constraints 
to  the  relative  velocities  at  the  entry  and  exit  points,  even  though  the  entry/exit 
/  are  free  parameters.  The  mission  scenario  might  be  for  an  observation  or  sentry 
satellite  that  is  already  in  a  near  orbit  with  its  period  matched  to  the  chief  but  needs 
to  “pause”  in  a  certain  quadrant  of  the  chief  in  order  to  perform  a  task  and  then 
continue  on  in  a  non-drifting  closed-relative  orbit  that  may  or  may  not  be  the  same 
as  the  one  at  entry.  Like  Case  1,  convergence  to  a  final  trajectory  happens  after  about 
two  or  three  legs,  therefore,  in  the  interest  of  brevity,  only  k  <  3  results  are  shown. 

This  particular  exit  condition  presents  a  harsher  scenario  for  the  continuous- 
thrust  solution.  At  the  start  of  the  exit  leg,  the  deputy  is  at  the  x^nm  position  with 
zero  relative  velocity.  The  algorithm  searches  for  a  minimum-fuel  /  from  which  to 
leave  the  lobe  under  the  condition  that  at  that  point,  the  deputy  must  thrust  in  order 
to  match  the  relative  velocity  of  a  closed  relative  orbit.  In  other  words,  the  deputy 
has  to  accelerate  from  zero  to  a  fairly  large  relative  velocity  that  may  occur  in  one  or 
two  separate  burns  (one  to  leave  and  one  to  exit  the  lobe).  On  the  other  hand, 
the  discrete  solution  does  not  slow  to  a  zero  relative  velocity  during  its  trajectory  and 
therefore  the  final  exiting  AV  is  smaller.  This  results  in  discrete-thrust  solutions  that 
have  significant  fuel  savings  over  the  continuous  solution  for  all  three  lobes. 

There  is  one  other  feature  of  note.  In  the  single  leg  of  Lobe  2  and  all  three 
examples  in  Lobe  3,  the  trajectory  seems  to  have  an  extra  leg.  This  is  a  result  of 
the  algorithm  finding  an  exiting  closed  relative  orbit  on  which  to  leave  that  continues 
to  pass  through  the  lobe  after  the  burn.  The  algorithm  detects  that  the  exit  leg 
remains  inside  the  lobe  and  adds  the  additional  time-of-flight  to  the  discrete  thrust 
solution.  For  example  in  Figure  55,  we  have  requested  a  single  leg  trajectory.  The 
algorithm  finds  '02  =  115°  as  the  optimal  exit  point  that  corresponds  to  a  closed 
relative  orbit  originating  on  the  upper  side  of  the  lobe  and  finally  exiting  on  the  lower 
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side.  The  calculated  time-of- flight  is  0.127  fractions  of  a  chief  period;  however,  the 
algorithm  determines  that  the  final  leg  remains  in  the  lobe  resulting  in  0.06  fractions 
of  additional  hover  time. 
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Figure  55:  Case  2,  Lobe  2,  Legs  =  1 
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Figure  56:  Case  2,  Lobe  2,  Legs  =  2 
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Figure  57: 


Case  2,  Lobe  2,  Legs  =  3 


(a)  XY  Plane  Trajectory 
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Figure  58: 


Case  2,  Lobe  3,  Legs  =  1 
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Figure  59: 


Case  2,  Lobe  3,  Legs  =  2 
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(a)  XY  Plane  Trajectory 
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Figure  60: 


Case  2,  Lobe  3,  Legs  =  3 
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Figure  61: 


Case  2,  Lobe  4,  Legs  =  1 


(a)  XY  Plane  Trajectory 
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Figure  62:  Case  2,  Lobe  4,  Legs  =  2 
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Figure  63:  Case  2,  Lobe  4,  Legs  =  3 
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5.4  Case  3:  Open  Entry  Condition/ Repeating  Hover  Orbit 

By  specifying  open  entry /repeating  hover,  we  seek  to  find  closed  trajectories 
that  can  be  sustained  for  as  long  as  the  mission  planner  desires  and  the  deputy’s 
fuel  stores  last.  This  is  done  by  constraining  the  relative  position  and  velocity  at  the 
final  point  (post  thrust)  to  be  the  same  as  the  initial  relative  position  and  velocity. 
These  closed  trajectories  are  recommended  when  simply  increasing  the  number  of  legs 
does  not  yield  the  needed  total  time-of- flight.  Two  types  of  closed-relative  orbits  are 
obtained  by  the  algorithm.  The  first  is  the  classic  teardrop  in  which  the  trajectory 
intersects  itself  at  a  single  point,  providing  an  opportunity  to  perform  an  impulsive 
thrust  and  repeat  the  teardrop.  The  second  is  called  the  bounce  trajectory  (Figure  65) 
and  occurs  when  the  deputy  bounces  back  and  forth  between  two  points.  The  teardrop 
trajectory  outperforms  the  bounce  and  is  the  focus  of  further  discussion. 
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Figure  64:  Case  3,  Lobe  2,  Legs  =  3,  Teardrop 
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(a)  XY  Plane  Trajectory  (b)  Growth  of  AV 


Figure  65:  Case  3,  Lobe  2,  Legs  =  2,  Bounce 
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Figure  66:  Case  3,  Lobe  3,  Legs  =  3,  Teardrop 
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Figure  67:  Case  3,  Lobe  4,  Legs  =  3,  Teardrop 
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Figure  68:  Case  3,  Lobe  4,  Legs  =  4,  Teardrop 
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Just  as  the  PHO  analysis  in  Section  5.1  provided  insight  into  a  more  efficient 
and  accurate  calculation  of  orbits  about  the  equilibrium  condition,  so  too  does  ex¬ 
amination  of  the  Case  3  trajectories  (Figures  64-68).  In  all  simulations,  the  lowest 
cost  solution  is  a  teardrop  with  expected  numerical  imprecisions  producing  orbits 
that  are  not  perfect  teardrops  returning  to  a  single  point  but  close  approximations. 
This  encourages  analysis  of  all  possible  teardrops  and  their  cycle  cost  {AV  per  cycle). 
Figures  69-71  present  the  results  of  that  analysis.  In  each  example,  only  ?/^’s  that 
have  T’max(t^i, t^i)  >  0  are  evaluated  (see  Figure  16,  Section  3.6  and  accompanying 
discussion)  and  only  for  T  <  T’max(t/'i,  t^i),  thus  producing  regions  of  no  data  in  the 
cost  surface.  That  cost  surface,  as  seen  in  sub-figures  (a),  show  that  although  there 
is  a  single  global  minimum,  a  family  of  nearly  identical  cost  teardrops  are  available 
to  a  mission  planner.  Thus,  if  cycle  period  is  an  important  mission  parameter,  a  wide 
range  of  values  can  be  selected  with  only  small  changes  to  cycle  cost.  Sub-figures  (b) 
display  a  representative  of  this  family  of  choices  chosen  for  its  larger  time-of-flight. 
Sub-figures  (c)  demonstrate  the  repeating  pattern  of  the  relative  specific  velocities. 

It  can  be  shown  (Appendix  N)  that  thrusting  at  the  apex  of  a  teardrop  in 
order  to  remain  on  that  relative  orbit  is  equivalent  to  thrusting  continuously  at  the 
time-averaged  x  coordinate  (centroid)  of  the  teardrop  (Equation  (144)). 


1  2(1-C)x„ 

^TD  =  75; —  /  x[t)dt  = 


T, 


TD  JO 


8  —  GtvTj'dS  —  SC 


(30) 


The  total  AV  per  cycle  can  be  found  closed  form  (Equation  (145)): 


AVtd  — 


127rrr£)(l  -  C)xo 
8  -  GttTtdS  -  SC 


(31) 


This  is  plotted  in  sub-figures  (d)  and  track  exactly  with  the  stairstep  of  the  discrete- 
thrust  solution.  Sub-figures  (d)  can  be  used  to  make  decisions  between  the  discrete 
and  continuous  solutions.  Note  that  these  solutions  assume  that  the  deputy  started  in 
the  desired  position  (start  of  the  teardrop  for  discrete  and  at  x^m  for  the  continuous) 
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Figure  69:  Lobe  2  Teardrop  Analysis 


with  the  desired  relative  velocity.  Fuel  usage  to  get  to  those  initial  conditions,  however, 
must  be  taken  into  account.  That  minimum- fuel  solution  is  treated  here  as  a  simple 
constant  added  to  a  linear  function  that  will  shift  the  AV  growth  profiles  up.  They 
will  not,  in  general,  be  the  same  for  the  continuous  and  discrete  solutions.  The  total 
fuel  used  in  the  discrete  solution  will  be 
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Tt 

Pd 


+  ^^D,c 


where  AVd  is  the  specific  AV  required  at  each  teardrop  burn,  Tt  is  the  total  time- 
of-flight  required,  Pd  is  the  period  of  the  teardrop,  and  AVdjc  required  to 
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Figure  70:  Lobe  3  Teardrop  Analysis 


get  into  position  for  the  first  teardrop.  The  continuous  solution  is 
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(a)  AV  Cost  Per  Cycle 
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Figure  71:  Lobe  4  Teardrop  Analysis 


This  is  our  criterion  for  choosing  between  a  discrete  and  continuous  solution.  In 
all  cases,  the  continuous  solution  will  eventually  outperform  the  discrete,  but  for  a 
shorter  total  time-of-flight,  it  makes  sense  to  use  a  discrete  strategy. 
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5.5  Trajectory  Changes  Due  to  Increased  xl 

As  noted  in  Section  3.6.7,  the  maximum  time-of-flight  between  any  two  de¬ 
creases  as  the  lobe  is  pushed  farther  from  the  Y  axis  (increase  in  the  lobe  center  x 
coordinate,  xl).  This  is  due  to  larger  relative  velocities  that  are  a  result  of  larger  dif¬ 
ferences  in  the  inertial  orbit  period  between  the  chief  and  deputy.  Figure  72  shows  a 
series  of  lobes  with  increasing  xl.  In  effort  to  compare  like  entry  and  exit  conditions, 
this  is  a  defined  entry/defined  exit,  both  of  which  have  relative  velocities  of  (0,0,0). 
Had  the  entry  velocity  been  chosen  based  on  a  closed  relative  orbit  (as  in  earlier 
examples),  the  velocity  would  have  changed  dramatically  as  increased.  Defining 
the  entry  velocity  as  (0,0,0)  provides  a  more  consistent  comparison  of  the  three  sim¬ 
ulations.  Two  legs  were  specified,  but  in  each  simulation,  the  second  burn  occurs 
nearly  simultaneously  with  the  exit  burn.  There  are  slight  changes  to  the  trajectory 
(both  discrete  and  continuous)  with  increasing  however,  more  importantly  is  the 
significant  change  in  time-of-fiight.  Again,  this  is  to  be  expected  after  examining  the 
effect  of  xl  on  the  constraint  surface  and  is  verified  here. 
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5.6  Trajectory  Changes  Due  to  Increased  Lobe  Size 

Changing  the  size  of  the  lobe  also  has  an  effect  on  trajectory  and  time-of- flight. 
Again,  from  our  discussion  in  Section  3.6.6,  we  know  that  increasing  the  overall  size  of 
the  lobe  increases  the  maximum  time-of- flight.  Figure  73  illustrates  the  sensitivity  to 
lobe  size.  The  changes  in  trajectory  are  more  pronounced  than  they  are  in  Section  5.5 
and  as  expectded  a  larger  time-of-flight  is  produced  with  larger  lobe  size. 
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Figure  73:  Defined  Entry/Defined  Exit,  Legs  =  2 
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5.7  Trajectory  Check  in  Other  Quadrants 

As  an  illustration  of  our  claim  that  lobes  in  the  first  quadrant  are  representative 
of  lobes  in  each  of  the  other  quadrants  (see  Section  3.6),  the  optimization  for  Case  1 
(Defined  Entry/Open  Exit)  with  two  legs  is  found  for  the  mirror  of  Lobe  2  (that  is, 
a  =  225°).  Figure  74  compares  the  two  trajectories  and  we  can  see  that  they  are,  as 
expected,  identical  once  the  rotations  discussed  in  Section  3.2.3  are  applied. 


Figure  74:  Comparison  of  Lobe  2  Results  to  its  Mirror  Lobe 


5.8  Addition  of  Z  Axis  Motion 

As  discussed  in  Section  4.2,  motion  in  the  Z  direction  is  completely  decoupled 
from  the  XY  plane  and  can  be  calculated  and  appended  to  the  XY  solution  inde¬ 
pendently.  A  few  examples  of  this,  based  on  previous  results,  are  provided  below. 
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Figure  75:  Case  1,  Lobe  3,  Legs  =  3,  /?  =  45°,  h  =  0.5  km 
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(a)  XY  Plane  (b)  3D  Trajectory 
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Figure  76:  Case  2,  Lobe  2,  Legs  =  2,  (3  =  45°,  h  =  0.5  km 


96 


(a)  xr  Plane 


(b)  3D  Trajectory 


Position  Histories 


Time  (fractions) 


Time  (fractions) 


(c)  Position  History 


Figure  77: 


Case  3,  Lobe  2,  Legs  =  4,  /3  =  45°,  h  =  0.5  km 
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5.9  Results  from  Elliptical  Chief  Orbit 

The  elliptical  chief  results  are  not  much  different  from  what  has  been  pre¬ 
sented  in  the  circular  case.  Recall  that,  unlike  the  circular  case  that  can  use  a 
single  constraint  surface,  the  elliptical  maximum  time-of-flight  is  a  function  of  ini¬ 
tial  true  anomaly.  This  means  that  practical  calculation  of  an  optimal  trajectory 
using  precomputed  constraint  surfaces  requires  a  library  of  those  surfaces  discretized 
over  true  anomaly.  The  following  result  is  for  a  single  leg  on  Lobe  2,  with  defined 
entry  {'ipi  =  45°)  and  open  exit  with  an  entry  velocity  of  (0,0,0).  It  demonstrates  that 
results  from  the  elliptical  case  can  be  found  with  no  change  to  the  cost  function  and 
are  similar  to  the  circular  case  (compare  to  Figure  40). 


Figure  78:  Case  1,  Lobe  2,  Legs  =  1,  Elliptical  Chief  Orbit 


5.10  Summary 

The  results  produced  above  show  that  there  are  a  wide  variety  of  practical 
lobes  available  to  mission  planners.  There  are  also  simple  tests  to  indicate  the  most 
fuel-optimal  method  (continuous  or  discrete)  to  use  for  a  specified  lobe,  entry/exit 
conditions,  and  required  total  time-of-flight.  It  is  also  important  to  note  that  there 
are  mission  scenarios  for  which  a  discrete-thrust  method  is  preferable  even  if  it  is  less 
fuel-optimal  than  the  continuous.  Issues  such  as  stealth,  plume  impingement  on  the 
chief,  and  sensor  vibrations  may  lead  a  mission  planner  to  the  discrete  strategy. 
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VI.  Conclusions  and  Recommendations 


The  goal  of  this  research  was  to  answer  two  questions.  The  first,  “Can  a  discrete 
trajectory  be  found  that  outperforms  the  benchmark  continuous-thrust  solution  for  a 
fuel  criterion  of  optimality?”  The  second  is,  “Can  we  quickly  and  robustly  estimate, 
with  reasonable  accuracy,  the  amount  of  AV  required  to  stay  in  a  specific  lobe?” 

The  answer  to  the  first  question  is  a  qualified  yes.  The  qualification  is  that  it 
can  be  accomplished  only  in  certain  combinations  of  entry /exit  conditions  and  then 
only  for  a  relatively  short  total  time-of- flight.  In  fact,  entry /exit  conditions  appear 
to  have  the  greatest  effect  on  which  strategy  is  fuel-optimal.  For  a  mission  planner 
attempting  to  decide  between  a  discrete-  and  continuous-thrust  solution  for  a  given 
lobe  and  entry/exit  conditions,  he  or  she  should  start  with  a  single  leg  and  find 
optimal  trajectories  for  about  k  <  5.  If  any  of  these  meet  the  total  time-of-flight 
required,  then  a  simple  check  of  the  AV  growth  charts  will  indicate  the  fuel-optimal 
method.  If  the  total  time-of-flight  cannot  be  satisfied  by  increasing  the  number  of 
legs,  an  evaluation  of  repeating  hover  orbits  is  needed,  as  discussed  in  Section  5.4,  in 
which  a  simple  test  can  determine  the  best  choice  of  longer-term  hovering  orbits.  As 
mentioned  in  Section  5.10,  there  may  be  other  concerns  besides  fuel-optimality  that 
may  lead  a  mission  planner  to  the  discrete  solution. 

The  answer  to  the  second  question  is  a  definite  yes.  The  continuous-thrust  so¬ 
lution  is  robustly  and  rapidly  calculated  without  any  pre-computing  of  the  constraint 
surface.  It  is  also  a  very  reasonable  estimate  of  the  total  AV  required  to  hover  in  a 
specific  lobe.  That  estimate  can  be  scaled  by  the  inverse  of  the  chief’s  mean  motion 
(n)  to  get  the  actual  AV  cost  for  a  given  sized  orbit.  The  ability  to  compare  all 
possible  size  chief  orbits  is  a  very  important  aspect  of  this  research. 

A  first  step  was  made  towards  applying  this  technique  to  relative  motion  about 
chiefs  in  elliptical  orbits.  Although  it  is  much  more  computationally  expensive  to 
calculate  the  relative  velocities  and  to  pre-compute  a  family  of  constraint  surfaces 
with  different  true  anomalies,  the  core  method  of  using  a  minimum- fuel  per  time-of- 
ffight  cost  function  and  a  maximum  time-of- flight  constraint  surface  is  sound. 
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Proximity  operations  and  hovering  specifically  are  fuel-expensive  operations  to 
perform  and  will  only  be  undertaken  when  mission  needs  are  great.  These  tasks  are 
not  however  out  of  the  AV  capability  range  of  current  micro-satellites. 

6. 1  Contributions 

Although  work  in  the  realm  of  relative  satellite  motion  and  proximity  operations 
has  been  extensive,  this  research  has  added  several  contributions  to  the  field.  Most 
significantly,  this  research  enables  mission  planners  to  search  for  and  compare  con¬ 
tinuous  and  discrete  fuel-optimal  trajectories  that  are  constrained  to  remain  within 
a  lobe  of  the  user’s  choosing  and  are  not  tied  to  a  particular  chief  semi-major  axis  or 
11.  In  completing  that  task,  several  minor  contributions  were  also  made: 

1.  Created  a  set  of  mission-realistic  lobe  entry  and  lobe  exit  conditions  for  use  in 
initializing  trajectories 

2.  Provided  an  analysis  of  maximum  time-of-flight  surfaces  and  their  sensitivity  to 
various  factors 

3.  Developed  a  simplified  optimization  routine  for  finding  persistent  hover  orbits 
that  intersect  the  Y  axis 

4.  Derived  a  closed-form  solution  for  optimal  hovering  along  the  Z  axis  that  is 
independent  of  the  trajectory  in  the  XY  plane 

5.  Applied  this  method  to  an  elliptical  chief  orbit  scenario 

6.  Provided  a  full  derivation  of  the  partial  derivatives  of  the  cost  function 

Mission  Impact:  This  research  enhances  USAF  proficiency  in  satellite  prox¬ 
imity  operations  and  space  situational  awareness  be  enabling  a  practical  method  of 
finding  relative  motion  trajectories  that  allow  hovering  near  a  chief  satellite. 
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6.2  Topics  for  Future  Research 

This  research  provides  only  the  first  step  into  hovering  orbits.  A  wide  array  of 
topics  are  available  to  future  researchers 

1.  Expand  the  work  on  hovering  about  chiefs  in  elliptical  orbits  to  include  newer 
research  on  more  efficient  relative  velocity  calculations. 

2.  Re-pose  the  cost  function  or  add  constraints  to  allow  the  possibility  of  thrust 
points  interior  to  the  lobe.  Compare  those  results  to  those  trajectories  produced 
when  thrusting  is  allowed  only  on  the  lobe  boundary. 

3.  Apply  additional  constraints  such  as  total  time-of- flight  or  minimum  AV  per 
burn  to  the  cost  function.  This  may  prevent  the  thrust  splitting  phenomena 
documented  in  Chapter  V  and  would  allow  inclusion  of  hardware  related  con¬ 
straints. 

4.  Include  attitude  control  in  addition  to  trajectory  planning  for  completeness  of 
the  total  mission  requirements.  What  if  we  wanted  the  deputy  pointing  towards 
the  chief  while  it  is  hovering? 

5.  Develop  a  strategy  to  include  additional  thrusts  that  correct  the  trajectory  for 
perturbations  and  equations  of  motion  linearization  errors.  This  is  especially 
important  for  repeating  hover  orbits  that  may  drift  out  of  the  lobe  over  time. 
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Appendix  A.  Math  Preliminaries 


This  appendix  contains  mathematical  theorems  and  proofs  that  are  used  in 
subsequent  derivations.  None  of  these  derivations  are  original  to  this  work  nor  partic¬ 
ularly  hard  to  hnd  in  mathematical  or  engineering  texts.  They  are,  however,  key  to 
later  derivations  and  included  here  for  completeness  and  to  put  them  into  dissertation 
notation. 


A.l  Theorem:  Radius  of  an  Ellipse 

The  radius  of  an  ellipse  measured  from  the  origin  is 

ab 

V sin  9  +  b'^  cos^  6 

with  all  variables  as  defined  in  Figure  79. 

Proof: 

Referencing  Figure  79,  choose  an  arbitrary  point  (x,  y)  on  an  ellipse  that  is  centered 


Figure  79:  Ellipse  Centered  on  the  Origin 
at  the  origin.  In  polar  form: 


X  =  r  cos  9 
y  =  r  sin  6* 
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Substituting  these  equations  into  the  standard  form  of  an  ellipse, 


solving  for  r, 


we  hnd  that, 


a?  62 

cos^  6  sin^  6 

- ^ —  + 


62 


=  1 


cos^d  sin^d"' 
— ^  + 


62 


=  1 


[6^  cos^  6*  +  sin^  6*]  =  o?b^ 
a%‘^ 


J.2  ^ 


62  cos2  6*  +  a2  sin^  0 


ab 


r  = 


\/62  cos2  d  +  a2  sin^  6 


(32) 


Theorem:  The  Inertial  Derivative 

The  derivative  of  a  vector  defined  in  a  rotating  reference  frame  is 

ijkdjy)  _  xYzd’jV)  \  ^ 
dt  dt 

where  to  is  the  angular  velocity  of  the  XY Z  frame  with  respect  to  the  IJK  frame 
and  the  notation  indicates  a  derivative  taken  in  the  ijK  frame. 

at 

Proof: 

The  following  is  a  modification  of  the  derivation  given  in  Wiesel  [63].  Let  IJK  be 
the  inertial  frame  and  XY Z  a  rotating  local  reference  frame  and  let  the  vector  V  be 
defined  in  the  XY Z  frame,  then 


V  =  KX  +  VyY  +  vj 


103 


Taking  the  derivative  using  the  chain  rule 


ijk^^ 


Noting  that 

=  14X  +  VyY  +  14Z 

The  changes  in  the  XY Z  frame  unit  vectors  can  be  found  by  referencing  Figure  80. 
Assume  the  local  frame  has  rotated  during  the  period  At.  This  rotation  is  formed  by 
inspection  of  Figure  80. 

X'  =  cos  (uAt)  X  +  sin  (uAt)  Y 
Y'  =  —  sin  (ojAt)  X  +  cos  (uAt)  Y 
Z'  =  Z 


Thus 


AX  =  X'  —  X  =  cos  {ujAt)  X  +  sin  (uAt)  Y  —  X  =  [cos  (cuAt)  —  1]X  +  sin  (cuAt)  Y 
AY  =  Y'  —  Y  —  —  sin  (cuAt)  X  +  cos  (uAt)  Y  —  Y  =  —  sin  (cuAt)  X  +  [cos  (uAt)  —  1]  T 

AZ  =  Z'  -Z  =  0 


As  At  goes  to  0 


uAt  —>■  udt 
AX  dX 

AY  dY 

cos  {(jjAt)  1 

sin  [uiAt]  ujdt 
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X 


then 


Figure  80:  Derivatives  of  the  XYZ  Frame  Unit  Vectors 


dX  =  XX  =  ujdtY 
dY  =  AY  =  -udtX 


Dividing  through  by  dt 


UK 


UK 


UK' 


dX 

dt 

dY 

dt 

dZ 


dt 


cuY 

-uX 

0 


Equation  (33)  becomes 

IJk^  =  Y-V,uiX  +  0Z 

dt  dt  ^ 

XYZ 

0  0  a; 

U.  Vy  K 


XYZ' 


dV 

dt 
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Therefore  the  inertial  derivative  is 


^  XYZ^y  ^  ^  X  y 
dt  dt 


(34) 


A. 3  Theorem:  The  Harmonie  Addition  Theorem 

The  following  derivation  is  taken  from  Reference  [61].  The  sum  of  two  sinusoids 
is  equal  to: 


a  cos  9  +  b  sin  6 
a  cos  6  +  b  sin  9 


+  6^ 
+  6^ 


cos 

sin 


e  -  ta.r‘  Q) 


Proof: 

Let 

f{9)  =  a  cos  9  +  b  sin  9 

Assume  we  want  to  express  this  as  a  function  of  a  single  sinusoid 


f{9)  =  ccos{9  +  uj) 
Expanding  the  above  equation  [53] 

f{9)  =  c  cos  9  cos  u  —  c  sin  6*  sin  a; 


which  means 


a  =  c  cos  UJ 
b  =  — csincu 


Finding  c 


2  ,  u2  2  2  ,  2  •  2  2 

a  +  b  =  c  COS  UJ  +  c  sm  uj  =  c 
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and  uj 


b  — csmcu 

-  = - =  —  tan  u 

a  c  cos  u! 


thus 


c  =  a?  +  h^ 


Therefore 


;  cos  6  -\-h  sin  6  =  +  6^  cos 


6  —  tan 


a 


We  can  also  express  this  in  terms  of  a  sine  function.  Let 


Expanding  [53] 


which  means 


f{0)  =  csin(6*  +  uj) 


f{9)  =  csindcoscu  +  c  cos  d  sin  a; 


(36) 


a  =  c  sin  cj 
b  =  c  cos  (jj 


Finding  c 


2,72  2-2  I  2  2  2 

a  +  b  =  csincL;  +  c  cos  u  =  c 


and  UJ 

a  c  sin  a; 

-  =  - =  tan  a; 

b  c  cos  UJ 


thus 


c 

UJ 


+  6^ 

tan-'  (5) 
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Therefore 


a  cos  9  +  b  sin  9  =  -v/a^  +  6^  si 


sin 


9  +  tan 


-1 


(36) 


A. 4  Theorem:  Partial  Derivative  of  u'RR'u  with  Respect  to  a  Scalar 

(u  is  a  Function  of  the  Scalar) 

The  following  derivation  is  introduced  in  [33] .  Let  i?  be  an  z  x  j  matrix  and  let 
the  components  of  u  be  functions  of  the  scalar  v;  then 

=  2n'RR'^ 

ov  ov 

Proof: 

Let 


A  =  u'RR'u  = 


Ui  U2  ...  Ui 


rii  ri2  ...  rij 

r2i  r22  ...  r2j 


Ui  ri2  ...  r 


1-3 


Til  r2i  ...  Til 

ri2  r22  ...  ri2 


Aj  r2j 


Ui 

U2 

Ui 


{uiTii  +  .  .  .  +  UiTii)  {uiri2  +  .  .  .  +  Uiri2)  .  .  .  {uiTij  +  .  .  .  +  UiVij) 


{uiTii  +  .  .  .  +  UiTii) 

{uiri2  +  . . .  +  Uiri2) 

{uiTij  +  .  .  .  +  UiTij) 


—  ('^1^11  +  .  .  .  +  UiTii)  +  {uiTi2  +  .  .  .  +  UiTi2)  +  .  .  .  +  {uiTij  +  .  .  .  +  UiTij) 
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Taking  the  partial  derivative: 


dA 

=  2  (miTii  +  +  •  • 


[  dui 

iTii)  rii 


+  2  (Miri2  +  U2r22  +  . . .  +  Uiri2)  (  ri2 


dv 

dui 

dv 


du2  , 
r2i^^  +  • 


^22' 


dv 

dU2 

dv 


dui 

■  +  v; 
dv 

dui 

dv 


=  2 


.  .  .  (  dux  du2  dui 

+  2  {uirij  +  «2r2i  +  . . .  +  UiTij)  (  ry—  +  rs,—  +  . . .  + 

Separating  into  matrix  form: 

(MiTii  +  .  .  .  +  UiT ii)  {uiri2  +  .  .  .  +  UiT  12)  .  .  .  {UiTij  +  .  .  .  +  U^T ij) 

and  separating  again 


=  2 


Ui  U2  ...  Ui 


Til  ri2 

^21  r22 

Til  ri2 


Aj 

r2j 

I'ij 


rii  r2i 
ri2  r22 

rij  r2j 


Til 

ra 

'f'ij 


dui 

dv 

du2 

dv 

duj 

_  dv  _ 

thus 


dlu'RR'u)  dA  ^  ,^^,du 

2'U  RR  77 

ov  ov  ov 


(37) 
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A. 5  Theorem:  Partial  Derivative  of  u'RR'u  with  Respect  to  a  Scalar 

(R  is  a  Function  of  the  Scalar) 

The  following  derivation  is  introduced  in  [33] .  Let  Rhe  an  i  x  j  matrix  and  its 
components  be  functions  of  the  scalar  v,  then 


d  {u'RR'u) 
dv 


2u'R 


dR' 

dv 


u 


Proof: 

Let 


rii 

ri2  ■ 

..  Tij 

rii 

ni  . 

•  •  rn 

Ui 

A  =  u'RR'u  =  ui  U2  . 

.  .  Ui 

^22  ■ 

..  r2j 

ri2 

1^22  • 

■  ■  ri2 

U2 

rn 

ri2  ■ 

..  Tij 

Aj 

ni  • 

..  Tij 

Ui 

(wirii  +  . . .  +  UiTii)  {uiri2  +  . . .  +  Uiri2)  . . .  {uRij  +  . . .  +  UiTjj) 


{uiTii  +  . . . 

{uiri2  +  . . . 

{uiTij  +  ... 


Taking  the  partial  derivative: 


dA 

=  2  {uiT II  +  U2r 21  + 


u. 


Til)  ^^1 


+  2  ('Uiri2  +  U2r22  +  . . .  +  UiVif)  (  Ui 


)^  +  . 

.  .  .  +  {uiTij  +  . 

.  .  T  UiTij 

dr2i 

dTii 

dv 

+  U2  ^  +.. 

dv 

■+Ui 

dv 

dri2 

dr22 

dTi2 

dv 

+  U2  ^  +.. 

dv 

■+Ui 

dv 

drij 

dv 

+  “2  o  +  •  • 

dv 

dv 

d-Uirn) 
+  Uiri2) 

+  UiTij) 


no 


Separating  into  matrix  form: 


=  2 


(tiirii  +  . . .  +  UiTii)  {uiri2  +  . . .  +  Uiri2) 


{uiTij  +  .  .  .  +  UiTij) 

(«1^ +  If) 


+  M  ^ 
■  ■  -  ^  dv 


and  separating  again 


=  2 


Vj\  112  •  •  • 


Til 

ri2  ■ 

..  rij 

dr  11 
dv 

dr  21 
dv 

drn 

dv 

Ui 

r2i 

r22  ■ 

■■  r2j 

dri2 

dv 

dr22 

dv 

dri2 

dv 

U2 

rn 

ri2  ■ 

..  rij_ 

dv 

dr2j 

dv 

drjj 
dv  _ 

Ui 

thus 


d  iu'RR'u) 
dv 


dA  ,dR' 

—  =  2u  R^u 
ov  ov 


(38) 
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Appendix  B.  The  User-Defined  Lobe 


Ideally  the  user  should  have  maximum  flexibility  in  specifying  the  lobe  in  which 
the  deputy  will  hover.  For  chiefs  in  circular  orbits,  the  out-of-plane  and  in-plane 
motion  decouple,  therefore  the  lobe  shapes  will  be  independently  constructed. 

B.l  In  the  Orbit  Plane 

In  the  orbit  plane  of  the  chief,  an  ellipse  is  a  convenient  closed  shape  that  pro¬ 
vides  the  utility  we  seek  without  overly  complicating  the  mathematics.  The  following 
is  a  derivation  of  the  polar  coordinates  for  an  arbitrarily  placed  and  oriented  ellipse. 
Reference  Figure  81  for  variable  deflnitions.  The  Cartesian  coordinates  of  the  deputy 


Figure  81:  The  Lobe  in  Two  Dimensions 
in  the  X'Y'  frame  are 


x'  =  r  cos  U 
y'  =  r  sin  U 


where,  from  the  proof  in  Appendix  A.l, 


r  = 


'^x'^y 


-y/ Ty  COS^  U  +  Tc  i’’ 


(39) 
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We  can  rotate  these  coordinates  into  the  chief  centered  frame  via 


X 

COS  77 

—  sin  77 

X' 

Y 

sin  77 

cos  77 

Y' 

thus 


Xi 

COS  77  —  sin  77 

Ti  COS 

r*  COS  77  COS  7/^'  —  Ti  sin  77  sin  'ip[ 

yi_ 

sin  77  cos  77 

Ti  sin  7/;' 

Ti  sin  77  cos  +  Ti  COS  77  sin  ip'^ 

Simplifying, 


Xi  =  Ti  cos  (77  +  'ifj'j) 
Vi  =  Ti  sin  [f]  +  'Ip[) 


Let  the  angular  position  {ipi)  be  measured  with  respect  to  a  line  parallel  to  the  X 
axis 

i’i  =  V  +  (40) 


Substituting 


Xi 


Vi 


y cos2  {ijji  -r])  +  T^  sin^  {ipi  -  rj) 

_ TxTy  cos  jjj _ 

y cos2  {tjji  -r])  +  T^  sin^  {ipi  -  rj) 

_ TxTy  sin  xpi _ 

^ r2  cos2  -ri)  +  T^  sin^  (7/;^  -  rj) 


(41) 


Translating  the  ellipse  so  that  it  is  centered  at  {xl^pl),  where 


xl  =  cos  «  =  7  cos  a  sin /3 
Pl  =  7^y  sin  a  =  7  sin  a  sin /3 
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yields 


.  ^  TxTy  cos  ibi  , 

Xi  =  'j  cos  a  sin  (3  H -  ^  (42) 

^ cos2  {'ipi  -r))+T^  sin^  {'ipi  -  77) 

.  ^  T^Ty  simpi  .. 

Hi  =  'J  sin  a  sin  (3  H -  ^  (43) 

J r2  cos2  {'ipi  -rj)  +  Tx  {A  -  v) 

if  Tx  =  Ty,  the  ellipse  becomes  a  circle  with  radius  r  and  the  equations  above  reduce 
to 


7  cos  a  sin  /3  +  r  cos  A 

(44) 

7  sin  a  sin  (3  +  t  sin  A 

(45) 

Next  we  need  the  partial  derivatives  of  Equations  (42)  and  (43)  with  respect  to  A- 
First  the  partial  derivative  of  the  denominator 


d  denom 

dA 


I  [-2Ty  cos{A  -  v)  sm{A  -v)  +  2tI  sm{A  -  v)  cos{A  -  v)] 

^ t2  cos2  {A  -V)  +  T^  sin^  {A  -  v) 

i'T'I  -  Ty)  cos(A  -  v)  sm{A  -  v)  _  n  {tI  -  T^)  cos{A  -  v)  sin(7l;j  -^^)^ 

t2  cos2  {A  -ri)  +  A  {A  -  v) 


Then 


Qx.  sj r2  cos2  {'(pi  -  77)  +  r2  sin^  {A  -  vAxTy{- sin  A)  -  '^xTy  cos  g” 

dA 

%  _  _ 

dipi  t2  cos2  (ipi  -v)+A  {A  -  v) 


cos^  {A  -  v)  +  A  {A  -  v) 

r2  cos2  {'ipi  -  77)  +  sin^  {A  -  vAxTyicos  A)  -  TxTy  sin  A^-^^ 
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Substituting  Equation  (41) 


dxj 

d^i 


X  y 


Sin 


i’i  -  (b?  -  cos(^j  -  if)  sin(^i  -  rf  cos  ipi 

X  y 


^x  — 

-Vi  sin  ^jji  —  rf  sin(2'0^  —  2r])  cos 


^  2r^r^ 

^  ‘  X  ‘  y 


(47) 


dVi 

dipi 


‘x  ‘y 


COS 


fi  -  Ti  ( 


cosifi  -  rf  sin(^j  -  rf  sin^j 


^2  2 
X  y 


Ti  COS  ipi  -  rf  sm{2ipi  -  2rf  sin  fi 

^^x'^y 


(48) 


and,  for  a  circular  lobe  (r*  =  =  Ty  =  r), 


dxj 

dipi 

dVi 

dfi 


—Tsinipi 
T  cos  ipi 


Conversion  back  to  polar  coordinates  from  Cartesian  is  setup  in  Figure  82.  Note  that, 


Figure  82:  Conversion  from  Cartesian  to  Polar  Coordinates 
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IxY  =  7sin/3 


The  relationship  of  x'  and  y'  to  xp  is 


tan^i  =  ^ 


XI 


Thus 


xpi  =  tan 


-1 


lA 

x'- 


=  tan 


-1 


l/j  —  7  sin  a  sin  f3 
Xj  —  j  cos  a  sin  (3 


and  the  radius  (r*)  is  given  by 


(49) 


f'i  =  V +  iVi^  =  \l {^i~l  COS  a  sin  (3f  +  {Vi  -  ^  sin  a  cos  (50) 


B.2  Out  of  the  Orbit  Plane 

In  three  dimensions,  the  lobe  is  constrained  to  be  an  elliptical  cylinder.  Em¬ 
ploying  a  one  dimensional  lobe  shape  in  the  out-of-plane  direction,  allows  us  to  take 
advantage  of  the  decoupling  of  the  equations  of  motion.  Thus  a  single  parameter,  (h) 
the  half  height  of  the  lobe,  is  required  to  dehne  the  lobe  in  three  dimensions. 


Figure  83:  The  Lobe  in  Three  Dimensions 
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Appendix  C.  The  Two  Body  Problem 


The  classical  two  body  orbit  problem  is  derived  below  for  completeness.  To  find 
the  acceleration  of  a  satellite  about  a  much  larger  primary  body  due  to  the  force  of 
gravity,  consider  Figure  84,  in  which  d  are  position  vectors  in  inertial  space: 


Figure  84:  Two  Body  Problem 
By  Newton’s  second  law  [38], 

F  =  Mass  *  Acceleration  =  ma  (51) 

In  the  presence  of  perturbations 


Fg  +  fc  +  fp  — 


ma 


(52) 


where  Fg  is  the  force  due  to  gravity,  p  is  a  vector  of  control  forces,  and  fp  is  a  vector 
of  perturbative  forces.  Gravity  obeys  an  inverse  square  law 


Fg\  OC 


mim2 
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where  ||c?||2  is  the  Euclidean  norm  of  the  relative  position.  To  make  this  proportionality 

/  .  ]XL^  \ 

equation  an  equality,  the  gravitational  constant  is  used  (  G  =  6.672  x  10“^^ - 5—  1 

V  kg  / 


I  p  I  _  Gmim2 

'  “  Ml 

In  order  to  apply  Newton’s  law,  we  need  to  move  the  inertial  reference  frame  so  that 
it  is  centered  on  the  larger  mass  (mi).  This  means  we  are  assuming  the  larger  mass  is 
not  accelerating  in  the  true  inertial  frame  and  thus  is  an  inertial  frame  itself.  As  long 
as  the  masses  are  sufficiently  disparate  and  there  are  no  other  forces  of  consequence 
operating  on  mi,  this  is  a  reasonable  approximation.  Substituting  into  Equation  (52) 
for  both  masses 


U  +  fci  +  fpi  =  midi 
{~u)  +  fc2-‘r  fp2  =  m2d2 

where  17  is  a  unit  vector  parallel  to  d  and  gives  the  scalar  force  of  gravity  a  direction. 
Making  the  following  substitution 


U  = 


d 


yields 


+  /cl  +  fpl  —  fTlidi 


+  fc2  +  fp2  —  m2d2 
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Combining  terms  and  dividing  through  by  the  masses  to  isolate  the  acceleration, 


di  = 


do  = 


Gm2 
—Gnii 


rui  mi 

m2  m2 


d  =  d2  —  di,  taking  the  derivative  twice  yields 


d  =  d2  —  di 


Substituting  Equations  (53)  and  (54)  into  the  above  equation 


—Gmi 

7  fc2  fp2  _ 

Gm2 

d-i^- 

_  Ui 

1 - 

toco 

1 _ 

(Jj 

m2  m2 

toco 

mi 

mi 

Combining  terms 


(53) 

(54) 


-G(mi  +m2)d  /d  fpi  ^  fc2  ^  fp2 


mi  mi  m2  m2 


If  the  first  body  is  the  Earth  and  the  second  body  a  satellite  orbiting  the  Earth,  then 

mi  dearth  5  ^2  ^sat?  fcl  /earth-controh  fc2  /sat-controh  fpl  /earth-pert  5  ^^id 
fp2  /sat-pert 


d  = 


+ 


G(7Tlearth  “1“  ^sat)^  /earth-control  /earth-pert 


\\d\ 


dearth 


dearth 


/e 


sat-control 

^sat 


/e 


sat-pert 

^sat 


Since  mearth  3>  TTisat  we  can  assume  that  mearth  +  ^sat  ^  dearth,  further,  if  the 
Earth  is  being  used  as  the  inertial  reference  frame  for  the  satellite  then  /earth-control  = 

/earth-pert  0. 


d  = 


G^^earthd  /sat-control  /sat-pert 


M 


^sat 


^sat 
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Defining  the  constant  ji  such  that 


/i 


N  • 

Gmearth  =  6.672  X  ^  X  5.974236  x  10^^  kg 


=  3.98601  X  10 


14 


m 


km 
=  398601- 


km^ 


and  defining  the  specific  forces  (accelerations) 


ttr  = 


CLj) 


/sat-( 


control 


^sat 

/sat-pert 


^sat 


Then  the  equation  of  motion  for  a  satellite  about  the  Earth  with  control  input  of  Uc 
and  a  perturbative  acceleration  of  is 


r  ~l^d 
cl  —  — z; —  “h  CLq  T 

Ml 


(55) 
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Appendix  D.  Orbital  Mechanics 


Several  orbital  mechanics  quantities  are  needed  in  subsequent  derivations.  This 
entire  appendix  originates  from  equations  and  discussion  found  in  [63]. 


Figure  85:  The  General  Elliptical  Orbit 


D.l  True  Anomaly  (v) 

Referencing  Figure  85,  the  PQW  is  an  inertial  frame  oriented  such  that  P  points 
towards  perigee,  Q  is  along  the  semi-latus  rectum,  and  IT  is  P  x  Q  along  the  angular 
momentum  vector.  The  position  vector  by  inspection  is 

r  =  r  cos  fP  +  r  sin  vQ 


Taking  the  derivative 

V  =[f  cos  V  —  rv  sin  P  +  [r  sin  z/  +  rh  cos  v]  Q 
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The  angular  momentum  vector  is 


H  =  fx  V  = 


Q 


w 


r  cos  u  r  sin  u  0 

r  cos  V  —  ru  sin  u  r  sin  p  +  ru  cos  u  0 

OP 

OQ 


(rr  cos  p  sin  p  +  r^z>  cos^  p) 
-{rf  cos  sin  1/  —  r‘^p  sir?  p)W 


=  r^pW 


The  angular  momentum  is  also  equal  to  [63] 


H  =  ^/jlp 


Substituting  r  from  the  next  section  (Equation  (64)) 


/ —  1  +  ecosi/  H?  r-,  ,  N 

rp  =  - =  Jjip  - =  \  ~  (1  +  ecosz^) 

r  p  y  p 


(56) 


Note  that  [63] 


Ip 

p 

p 


rp 


I  +  e  cos  p 
=  ^iSMA(l  —  e^) 


n  = 


p 


‘SMA 


(57) 

(58) 

(59) 


Solving  for  p  with  another  substitution  of  r, 


IP  ,1  +  ecosz/ 

p  =  J—  (1  +  ecosz^) - 

p  p 

~p  (1  +  ecosi^)^ 


(1  +  ecosz^)" 


o; 


SMA 


(l-e2)i 


n(l  +  ecosi^)^ 
(l-e2)f 


(60) 
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The  angular  rate  can  also  be  expressed  as 


r^z>  = 

.2  l-cp  fir(l  +  ecosu)  fj,(l  +  ecosu) 


Taking  the  derivative  of  Equation  (60) 


h  =  - ^2(1  +  ecosz/)e(— sinz/)h 

(1  -  e2)2 

— 2nez>(l  +  ecosz/)  sinz/  — 2ez>^sinz/ 
_  g2^|  1  +  ecosu 


(61) 


(62) 


or,  in  terms  of  only  z/, 


z> 


— 2esinz/  n(l  +  ecosz/)^ 
1  +  ecos  u  (1  —  e2)f 

— 2n^e(l  +  e  cos  z/)^  sin  u 

(1  -  e2)3 


n  2 


(63) 


£),;2  Ellipse  Radius  (r) 

The  radius  of  an  orbit  is  a  function  of  true  anomaly  (z^)  and  comes  from  the 
equation  of  a  conic  section  [63] 


Qsma(1  -  e^)  _  P 
1  +  e  cos  1  +  e  cos  u 


(64) 


where  p  is  the  semi-latus  rectum.  The  derivative  is 


r  =  p  [— (1  +  ecosi^)  ^(— z>esinz/)]  = 


peiy  sm  z/ 


e  sm 


=  rp- 


(1  + ecos  1^)2  1  +  ecosi^ 


(65) 


Substituting  Equation  (56) 

[P  . 

T  =  ^  —  e  sm  p 

V  p 


(66) 
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and  finally  the  second  derivative  of  r 


r  = 


ev  cos  V 


and  substituting  Equation  (57) 


r  = 


ru 


1  +  e  cos  u 


ev  cos  V 


cos  v 
1  +  e  cos  V 


(67) 
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Appendix  E.  General  Relative  Motion  Derivation 

Assume  the  chief  satellite  is  in  a  closed  Keplerian  orbit  and  define  the  following 
chief  centered  reference  frame:  the  X  is  oriented  along  a  line  from  the  center  of  the 
Earth  to  the  chief,  Z  is  perpendicular  to  the  orbit  plane  of  the  chief  and  Y  completes 
the  frame  as  the  cross  product  Z  x  X.  The  in-track  direction  is  aligned  with  the 
velocity  vector  of  the  chief  when  in  a  circular  orbit.  This  frame  is  commonly  referred 
to  as  the  Local- Vertical/Local-Horizon  (LVLH)  frame.  By  inspection  of  Figure  86, 
the  position  vectors  in  the  LVLH  frame  are: 

L  =  roX  +  0Y  +  OZ 
P  =  xX  -\-  yY  +  zZ 

M  =  L  +  P  =  {to  +  x)X  -[-  yY  +  zZ  (68) 

where  is  the  instantaneous  orbit  radius  of  the  chief  satellite  and  is  a  function  of 
time. 


(a)  Inertial  Vector  Relationships 
Figure  86: 


(b)  The  Chief  Centered  Relative  Frame 
General  Relative  Motion  Setup 


The  inertial  derivative  of  a  vector  written  in  a  rotating  frame  is  derived  in 
Appendix  A. 2  as 


V  = 


XYzd(X) 

dt 


-\-  oj  X  V 
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where  cu  is  the  angular  velocity  of  the  XY Z  frame  with  respect  to  the  inertial  frame. 
For  the  orbit  problem,  u  =  OX  +  OY  +  uZ  where  z>  is  the  time  rate  of  change  of 
the  chief’s  true  anomaly.  Taking  the  inertial  derivative  of  Equation  (68)  yields  the 
relative  velocity: 

M  =  +  {pZ)  X  M 

where  is  the  derivative  of  M  in  the  XY Z  frame. 


={fo  +  x)X  +  yY  +  zZ 


A 

Y 

z 

-pyX 

0 

0 

z> 

= 

p{ro  +  x)Y 

To  +  X 

y 

OZ 

Adding  these  together  yields: 


M  =  {to  +  X  —  py)X  -\-[y  +  ^{ro  +  x)]Y  +  zZ  (69) 

Taking  the  inertial  derivative  of  Equation  (69)  yields  the  acceleration: 


M  =  +  (uZ)  X  M 


where  is  the  derivative  of  M  in  the  XY Z  frame. 

=[r„  +  x  -  {yy  +  vy)]  X  +  [y  +  v{ro  +  x)  +  v{fo  +  h)]  F  +  zZ 


X 

Y 

z 

-  \py  +  p^{ro  +  x)]  X 

0 

0 

V 

= 

[pfo  -\-  px  —  p'^y]  Y 

fo  +  x-py 

y  +  p{ro  +  x) 

z 

- 1 

o 

Adding  these  together  yields: 

M  =  [fo  +  x  —  2py  —  i)y  —  p‘^{ro  +  a:)]  X+  [y  +  2h(ro  +  h)  —  +  u{ro  +  x)\  Y  +  zZ 
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Note  that  (via  substitution  of  Equations  (62)  and  (65)) 


2i'fo  +  PTo  =  2h 


1  +  e  cos  V 


2ez>^  sin  u 


1  +  e  cos  p 


fo  =  0 


and  also  (via  substitution  of  Equations  (61)  and  (67)) 


To  - 


1  +  e  cos  p 

jiil  +  e  cos  p) 


-  p^To  =  P^To 


1  +  e  cos  p 


e  cos  p 
1  +  e  cos  p 

—  =  ZR 


-  1 


-P^Tn 


1  +  e  cos  p 


Thus 


M  = 


r.  •  •  ••  -9 

X  —  zpy  —  py  —  p  X - - 


X  +  [y  -\-  2px  —  p^y  +  px]  Y  +  zZ  (70) 


We  know  from  Appendix  C  that  acceleration  due  to  gravity  and  control  forces  is: 


M  = 


—jiM 

PI 


+  Uc  T  Q-p 


where  y  is  the  gravitational  constant,  Ec  is  the  control  acceleration,  and  is  the 
perturbative  acceleration.  Expanding  the  magnitude  of  the  position  vector  yields 


IIMII2  =  a/ {to  +  xy  -\-y‘^  +  =[{ro-\-xf +  y‘^  +  z^] 


Thus 


M  = 


-lA{^o  +  x)X  +  yY  +  zZ] 
[{to  +  xy  +  y‘^  +  z'^]^ 

Setting  Equations  (70)  and  (71)  equal  to  each  other. 


■j-  Clr  Clj, 


(71) 


^ .  ..  .9 

X  —  2vy  —  vy  —  v  x - - 

X 

{To  +  x)X 

^c»_ 

 -y 

[y  +  2z>x  —  if^y  +  vx\  Y 

[{To  +  xY  +  +  z‘^]i 

yY 

+ 

zZ 

zZ 

L 

i^Cx  T  ^Px)^ 
i^Cy  +  ^Py)^ 
{ttcz  +  0,pz)^ 
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In  scalar  form  and  simplifying 


X  —  2i'y  —  i)y  —  —  [i 

y  +  2i'x  i)x  —  u'^y  + 


1 


X 


[{'f'o + + -2^]  2 


^Cx 


^Px  0 


[{to  +  xy  + 

fiz 


3  ^Cy  O'Py  0 


[{To  +  x)2  +  1/2  +  z2]2 


3  OjCz  ^Pz  0 


From  Equation  (61)  we  can  set  /i  equal  to 


fi  = 


r3jy2 


1  +  e  cos  u 


Substituting 


X  —  2uy  —  i)y  —  v'^x  — 


1 


y  +  2vx  -\-  i)x  —  u^y  + 


1  +  ecosuJ  [r2  +  x)^  +  y^  +  z^]\ 

;,2.,  ^  I  ^  _ y 


3  ^Cx  ^Px  0 


z  + 


1  +  e  cosu  J  _|_  _|_  ^2  _|_  ^2j  2 

d 


^Cy  ^Py  ~  0 


l  +  ecosz/y  [(^^  +  ^)2  +  ^2  +  ^2]2 


3  (^Cz  ^Pz  0 


Factoring  out  a  z>2  yields  the  general  relative  equations  of  motion 


X  —  2i)y  —  i)y  —  u 


;.2 


X  + 


To 


rl{ro  +  x) 


1  +  e  cos  V 


(1  +  e  cos  u)  [{Vo  +  x)2  +  1/2  -|-  z'^]  2 


3  ^Cx  ^Px  0 


(72a) 


y  +  2i'x  +  i)x  —  u^y 


1  - 


2  +  u  z 


(1  +  e  cos  u)  [{to  +  xY  +  ?/2  +  z'^]  2 
^3 


-  a.cy  -  «Py  =  0  (72b) 


(1  +  e  cos  v)  [{To  +  xY  +  ?/2  +  2 


^Cz  ^Pz  0 


(72c) 


We  can  further  simplify  these  equations  by  assuming  that  the  deputy  is  close  to  the 
chief  compared  to  the  instantaneous  orbit  radius  (ro).  Define  the  nonlinear  term  (NL) 
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as 


NL  =  [{vo  +  a:)^  + 


Expanding  and  factoring  out  an  r, 


NL  = 


2x  y 


2  ^2 


IH - 1 — ^ 


N  rl  rl 


2x  x"^ 


v2\  I 


—  r„(  IH - 1 — -  -\ — -  -\ — - 

To  rl 


Using  the  binomial  series  expansion  [24],  this  becomes 


NL  =  rt 


3  f  2x  x^  y‘^  1 


Neglecting  terms  of  order  larger  than  two,  0{x,y,z),  and  assuming  that  x,  y,  and  z 
are  appropriately  small  compared  to  the  instantaneous  radius  of  the  chief’s  orbit  (our 
underlying  assumption  for  linearization)  such  that 


y'^ 

_  (l 

^  ^  ^  U 

rpZ  rpZ  rpZ 

O  O  O 


the  nonlinear  term  is  now 


o  \  '  /  o  . 

To  J  \  To 


and  the  general  relative  equations  of  motion  are 


X  —  2uy  —  i)y  — 


x 


rl{ro  +  x) 


1  +  e  cos  u  (1  +  e  cos  v) 


^Px  0 


y  +  2ux  +  i)x  —  p'^y 


z  +  P^z 


1  - 


(1  +  e  cos  p)  rl  j  j 


^Cy  %y  ~ 


(1  +  e  cos  p)  rl 


^Cz  ^Pz  0 
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Canceling 


X  —  2uy  —  i)y  — 


X 


ro{ro  +  x) 


y  +  2ux  -\-  i)x  —  u'^y 


1  +  ecosi/  (1  +  ecos i/)  (fo  +  3a:) 

Tn 


^Px  ~  0 


1  - 


2  +  u^z 


(1  +  e  cos  p)  (fo  +  3a:) 

Tn 


Vy 


(1  +  e  cos  p)  (vo  +  3x) 


^Cz  ^Pz  0 


Combining  terms 


X  —  2py  —  py  —  p"^ 
y  +  2px  +  px  —  p^y 
z  +  Pz 


x  + 

1  - 


2xrn 


(1  +  e  cos  p)  (to  +  3x) 

_ To _ 

(1  +  e  cos  p)  {to  +  3a:) 

Tn 


^Px 


O'Cy  CCpy 


(1  +  e  cos  p)  (to  +  3x) 


^Cz  0,pz  0 


Assume  that  Vo  is  appropriately  large  compared  to  x  such  that 


=  0 
=  0 


To  +  3a:  ~  To 


Then 


X  —  2py  —  py  —  p^ 
y  +  2px  -\-  px  —  p‘^y 
z  +  p'^z 


x  + 

1  - 
1 


2x 


1  +  e  cos  p 
1 

1  +  e  cos  p 


^Cx  ^Px 


O-Cy  Ctpy 


1  +  e  cos  p 


^Cz  ^Pz  0 


0 

0 
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and  finally,  combining  terms,  the  simplified  general  relative  equations  of  motion  are 


X  —  2i'y  —  i)y  —  u‘^x 
y  +  2ux  +  ux  —  u'^y 


3  +  e  cos  V 
1  +  e  cos  u 
ecosv 
1  +  e  cos  u 


z  +  u‘^z 


ecosp 


-  a. 


Cz 


^Px  ~  0 

(73a) 

“  ^Py  ~  0 

(73b) 

=  0 

(73c) 

where 

n(l  +  ecosz/)^ 

^  =  - 3 — 

(l-e2)i 

— 2eh^sini/  — 2n^e(l  +  ecosz^)^sini/ 

1  +  e  cos  u  (1  — 

If  the  chief  satellite  is  in  a  circular  orbit,  then 


e  =  0 

if  =  n 

i>  =  0 


and  the  general  relative  equations  of  motion  (Equation  (72))  reduce  to 


X  —  2ny  —  v?{x  +  r^,) 
y  +  2nx  —  v?y 


1  - 


[{To  +  x)2  +  ?/2  _|_  22]2  ^ 
^3 


^Cx  ^Px  ~  0 


1  - 


5  +  n^z 


[{to  +  x)2  +  y‘^  + 

^3 


^Cy  %y  ~ 


_[(ro  +  x)2  +  1/2  -I-  2;2]  2  ^ 


^Cz  ^Pz  0 


(74a) 

(74b) 

(74c) 
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The  linear  form  of  Equation  (73)  reduce  to  the  familiar  Clohessy-Wiltshire  equations 


2ny  -  3n‘^x  -  =  0 

(75a) 

y  T  2nx  0,^^  ^py  ® 

(75b) 

z  +  n^z  —  ttc^  —  o-pz  =  0 

(75c) 
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Appendix  F.  Chief  Orbit  Fractions  as  the  Independent  Variable 

There  are  advantages  to  using  chief  orbit  fractions  as  the  independent  variable 
as  opposed  to  absolute  time,  primarily  the  disentanglement  of  the  relative  trajectory 
to  semi- major  axis  (or  orbital  radius  for  e  =  0).  The  relationship  between  absolute 
time  (t)  and  time  as  a  fraction  of  orbit  period  (i)  is 


where  the  chief  orbit  period  (P)  is  [63] 


P  =  27: 


27r 

n 


therefore 


n 


and  the  derivative  of  t  with  respect  to  t  is 

dt  2tt 
di  n 

The  relative  positions  remained  unchanged  in  this  conversion 


(76) 


X  =  X 

(77a) 

y  =  y 

(77b) 

z  =  z 

(77c) 

The  relative  velocities  as  functions  of  chief  orbit  fractions  are 


dx 

dx 

dx  dt 

2n 

(78a) 

X  = 

di 

di 

dt  di 

=  X  — 

n 

dy 

_  dy 

dy  dt 

2tt 

(78b) 

y  = 

di 

di 

dt  di 

=  y— 

n 

dz 

dz 

dz  dt 

.  Stt 

(78c) 

5  = 

di 

di 

dt  di 

n 
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and  the  accelerations 


dx  27r  dx 

di  n  di 

dy  2%  dy 

di  n  di 

dz  27r  dz 

di  n  di 


2w  dx  dt  .. 
n  dt  di  n? 

2tt  dy  dt 
n  dt  di  ^  'n? 

2tt  dz  dt  ..dTT^ 

n  dt  di  ri^ 


(79a) 

(79b) 

(79c) 


Also  let 


p  =  p 

^  .271  27r(l  +  e  cos  z/)^ 

p  =  p —  =  — ^ 

n  (l-e2)f 

”  ..dyr^  — 87r^e(l  +  e  cos  i/)^  sin 

^  ^  (1  -  e2)3 


We  are  now  ready  to  formulate  the  general  relative  equations  of  motion  as  functions 
of  chief  orbit  fractions.  Starting  with  the  homogeneous  form  of  Equation  (73)  from 
Appendix  E  and  dividing  through  by 


X  py  i)  p^ 
^  -  2— - -y - -X 


y  px  p  p 
^  +  2^  H — 5-x - -y 


z  z>^ 

^2  ^2^ 


1 


1  +  e  cos  p 


3  +  e  cos  p 
1  +  e  cos  p 
ecos  p 
1  +  e  cos  p 

=  0 


0 

0 


Substituting  the  new  acceleration,  velocities,  positions,  and  derivatives  of  true  anomaly 


1  ■■  n?  2  ■  n  ■  n  1  ■■  'n?  1  /  •  n  \  2 

'n?  dvr^  'n?  27r  27r  dzr^  V  27r/ 


3  +  e  cos  p 
1  +  e  cos  p 


- V 

n?  dvr^  n?  2tv  27r  dvr^  \  27r/ 


e  cos  p 
1  +  e  cos  p 


1  -  n? 

np.  47^2 


1 

1  +  e  cos  p 


=  0 


0 


0 


13d 


Simplifying  and  multiplying  through  by 


=  0  (80a) 

=  0  (80b) 

(80c) 

When  the  chief  is  in  a  circular  orbit 

e  =  0 
9  =  2tt 

i)  =  0 


X  —  2uy  —  uy  —  i)‘^x 
y  +  2i)x  +  ux  —  9‘^y 
z  +  i>‘^z 


1  +  e  cos  9 


3  +  e  cos  9 
1  +  e  cos  9 
ecos9 
1  +  e  cos  9 

=  0 


thus  the  Clohessy-Wiltshire  equations  a  functions  of  chief  orbit  fractions  are 


X  —  Ai^y  —  12tt‘^x  =  0 

(81a) 

y  +  Attx  =  0 

(81b) 

5  +  47r^5  =  0 

(81c) 
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Appendix  G.  A  Closed-Form  Solution  of  the  Linear 
Clohessy-  Wiltshire  Equations 

Start  with  the  homogeneous  linear  Clohessy-Wiltshire  equations  (Equation  (75) 
from  Appendix  E) 


X  —  2ny  —  3rSx  =  0 
y  +  2nx  =  0 
z  +  Fz  =  0 


Taking  the  Laplace  transform 

[s^A(s)  —  sXo  —  xf\  —  2n[sY(s)  —  yo]  —  3FX{s)  =  0 
[sW (s)  -  syo  -  ijo]  +  2n[sX{s)  -  =  0 

[s^Z(s)  —  szo  —  io]  +  n^Z{s)  =  0 

where  Xo  and  x^  are  initial  conditions  of  the  relative  position  and  velocity  in  the  X 
(radial)  direction,  yo  and  ijo  are  in  the  Y  (velocity)  direction,  and  Zo  and  Zo  are  in 
the  Z  (out-of-plane)  direction.  Collecting  terms  and  placing  the  equations  in  matrix 
form 


—  SF  —2ns  0 

.Y(s) 

SXo  ^Xo-  2nyo 

2ns  s^  0 

F(s) 

= 

syo  +  yo-  ‘^nxo 

0  0 

Z(s)_ 

SZq  Zq 

Einding  the  inverse  and  solving  for  X(s),  T(s),  and  Z{s) 


1 

2n 

X(s) 

-|-  n^ 

s{s^  -|-  n^) 

(J 

SXo  YXo-  2nyo 

—2n 

s^  —  SF 

Y(s) 

— 

s(s2  -|-  n^) 

n^) 

0 

syo  +  yo-  ‘inxo 

Z(s) 

0 

0 

1 

SZq  Zq 

136 


Multiplying  this  out 


sXo  +  Xo  —  2nyc 

,  2n{syo  +  ijo  + ‘2nXo) 

X(s) 

r(s) 

/<*>. 

+  n? 

—2n{sXo  +  Xo  —  2nyo) 

s{s‘^  +  ’n?) 

{s^ -Zv?){syo^-yo  +  2nXo) 

s{s‘^  +  n^) 

s^{s^  +  ’n?) 

SZq  T  Zq 

-\-  in? 

Performing  partial  fraction  expansion 


X(s) 

Y(,s) 

= 

Z(s) 

sXq  +  Xo-  2nyo  ^  2n‘^yo  -  2syo  -  AsnXp  ^  2yo  +  AnXp  1 
+  ri^  n{s^  +  n^)  n  s 

4sr/o  +  4yo  +  8nXo  —2n^Xo  +  2sXo  —  '^syo'n  Si/o  —  6nXo  nyo  —  2xo  1 


+  ri^ 


+ 


n(s2  +  V?) 

SZq  ~\~  Zq 


+ 


+ 


n 


+  'n? 


Collecting  terms 


X(s) 

y(s) 

Z{s) 


iih  , 

I  o  o  ' 

n  s 


+  ri^ 


2(^)s  2  (2yo±^  _  ^  2n  2%  +  4nXo  1  nyo  -  2xo  1 


+ 


+  n? 


SZo 


+ 


2  n 


+ 


n 


Taking  the  inverse  Laplace  transform 

2yo  +  4nXo 


x{t) 

y{t) 

Xr 


Xn  - 


2  (  ^  )  cos(nt)  +  2 


n 

2yo  +  4nx, 


.  ,x  *0  2yo  +  ^nXo 

cosint)  H - sinlnr)  H - 

n  n 


Xn 


n 


.  3nt  /2yo  +  4na:o\  nyo-2xo 


;2o  cos{nt)  H — -  sin(nt) 


n 


Since  the  ICs  and  n  do  not  vary  with  time,  two  constants  can  be  dehned 

2yo  +  ^nxo 


a  = 


6- 


n 

nyo  -  2xo 


n 


(82) 

(83) 

(84) 
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which  simplifies  Equation  (82)  to 


X 

x(t)  =  (xo  —  a)  cos(nt)  H — -  sin(nt)  +  a 

n 

t7i  n 

y{t)  =  — -  cos{nt)  +  2(a  —  Xg)  sm{nt) - —t  +  b 

Til  Zi 

zit)  =  Zg  cos(nt)  +  —  sinint) 

n 


(85a) 

(85b) 

(85c) 


For  simplicity,  x  =  x{t),  y  =  y(t),  z  =  z{t).  Isolating  the  sinusodal  terms. 


Xr 


x  —  a  =  (xg  —  a)  cos( nt)  H — -  sin(nt) 

n 

y  +  -b  Xg  ,  .  ,  ^  ^ 

- 7- - =  —  cos[nt)  +  [a  —  Xg)  sm{nt) 

Z  Th 


Squaring  both  sides  yields: 


[x  —  of  =  (xo  —  a)^  cos^(nt)  +  ^— J  sin^(nt) 

+  2(xg  —  a)  —  cos(nt)  sin(nt) 
n 


[y  +  ^t-  bf 


Xr 


n 


cos^  (nt)  +  (a  —  Xgf  sim(nt) 


—  2(xg  —  a)  —  cos(nt)  sin(nt) 
n 


Adding  these  equations  together 


Jy+lft-bf 

\x  —  a\  H - A - 


(x„  -af+{^ 

'  n 


cos^  (nt) 


+ 


+(-l)^(x„-a)^ 
n 


sin^(nt) 
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Pulling  out  the  common  coefficient  and  using  the  identity  cos^(q;)  +  sin^(Q;)  =  1 


X 


12  ,  [y 

o  H - 


+  ^t-  bf 


{xo  -  af  + 


2 


The  right  side  of  this  equation  has  no  time  varying  terms,  thus  define  the  constant  p 


=  (x,  -  af  + 

Substituting  and  dividing  by  yields 

[x  -  af  [y+^t-  hf 

p2  4p2 


(86) 


(87) 


This  is  the  standard  equation  for  an  ellipse,  thus  the  deputy  satellite  will  follow  an 
elliptical  path  about  the  chief  satellite  in  the  XY  plane  (Figure  87).  Note  that  for 
a  stable  orbit  (i.e.  one  that  does  not  drift  over  time)  about  a  chief  satellite  a  must 
equal  zero. 


Figure  87:  Relative  Orbit  Types 
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We  can  further  simplify  Equation  (85)  using  harmonic  addition  theorem  derived 
in  Appendix  A. 3. 


x{t)  =  {xo  —  a)^  +  (  —  )  sin 


nt  +  tan 


f  n{xo  —  a) 


Xo 


+  a 


y{t)  =  \l  (2)2  (  —  )  +  (-2)2  {xo  -  a)  cos 


nt  +  tan 


f  —2n(a  —  Xo) 

2Xn 


3na 


t  +  b 


Substituting  Equation  (86)  and  dehning  a  new  constant  0 


x{t)  =  psin  {nt  +  9)  +  a 


(88) 


^TiQj 

y{t)  =  2pcos  {nt  +  9) - —t  +  b 


where 


9  =  tan 


-1 


n{xc 


a 


Xn 


(89) 


(90) 


If  we  visualize  the  elliptical  orbit  from  Figure  87(a)  rotating  about  the  major  or  minor 
axes,  it  is  easy  to  see  it  will  trace  a  line  m  XZ  plane,  if  rotated  about  the  major  axis, 
and  the  Y Z  plane,  if  rotated  about  the  minor  axis.  Using  the  equation  of  a  line  and 
adding  the  effects  of  a  rotation  about  both  major  and  minor  axes 


z{t)  =  l[x{t)  —  a]  +  g 


,  ,  3na 

y{t)  +  ^t-b 


where  I  and  q  are  the  slopes  of  the  lines  formed  by  the  rotation  about  the  minor 
and  major  axes,  respectively.  Substituting  Equations  (88)  and  (89)  into  the  above 
equation  yields: 

z{t)  =  lpsm{nt  +  9)  +  2qpcos{nt  +  9)  (91) 

This  is  the  correct  form  expected  from  the  inverse  Laplace  transform  found  in  Equa¬ 
tion  (85c),  an  oscillatory  function  with  a  constant  amplitude  and  a  period  of  ^  or 

/  3 

2'kJ'^,  which  is  the  same  as  the  orbital  period  of  the  chief  satellite.  To  find  expres- 
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sions  of  I  and  q  in  terms  of  initial  conditions,  we  start  with  Equation  (88)  and  its 
derivative  (both  at  t  =  0) 


Xo  —  a  =  psm(9) 

Xo  =  pncos{6)  ^  ^  =  pcos(d) 

and  inserting  them  into  Equation  (95)  and  its  derivative  (also  at  t  =  0) 

Zo  =  lpsm(9)  +  2qpcos(9)  =  l(xo  —  a)  +  q  {  2— 

\  n 

Zo  =  lpncos{9)  —  2qpnsm{9)  =  Ixo  —  q  [2n{xo  —  a)] 

In  matrix  form 


~  ~ 

^0 

Xo  —  a  2 — 

1 

= 

n 

Zo 

Xo  —2n{xo  —  a) 

Q 

Solving  for  I  and  q  by  multiplying  both  sides  by  the  matrix  inverse 

—n^(a—Xo) 


Xq 


x‘^—2xon‘^  a+n'^  Xx"^  71“^  x‘^—2rn‘^  Xx‘^ 


^XoTl 


^{a-Xo)n 


_n‘^ x‘^—2xon‘^ a-\-n‘^ a‘^ -\-x^  n‘^x‘^—2xon‘^a-\-n‘^a‘^-\-x‘^  _ 

Multiplying  through  and  simplifying 

ZgXo  —  Zo'n?{a  —  Xo) 


I  = 


Q  = 


x^  +  n‘^(a  —  XoX 

ZoXoTi  +  Zon{a  —  Xo) 

2[xl  +  n'^(a  —  Xo)^] 


(92) 

(93) 


The  Z  motion  can  also  be  formulated  in  the  alternative  way.  Starting  with  Equation 
(85c) 

z{t)  =  Zo  cos(nt)  H — -  sm{nt) 
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the  sinusoidal  terms  can  be  combined  via  the  harmonic  addition  theorem  (Appendix  A. 3) 


z{t) 


+  Zg  cos 


nt  +  tan  ^ 


(94) 


Define  two  new  constants 


(f)  = 


2 


+  -2 


2 

o 


nZn 


Thus 


«(*)  =  Jii,.xCos(n(  +  tf)) 


(95) 


Using  the  conversions  formulated  in  Appendix  F,  we  can  express  these  closed-form 
solutions  as  functions  of  chief  orbit  fractions  as  opposed  to  absolute  time.  Noting  that 
nt  =  27rt  (Equation  (2)  in  Section  3.1) 


x{i)  =  x{i)  =  psin(27rt  -\-  6)  ^  a 
y{i)  =  y{t)  =  2pcos(27rt  -\-  9)  —  Sirai  +  b 
z{i)  =  z{i)  =  /psin(27rt  -|-  6*)  -|-  2gpcos(27rt  -|-  0) 


or 


Z{i)  =  5max  COs(27rt  -|-  0) 


(97) 
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where  the  values  of  the  relative  orbit  elements  do  not  change  when  calculated  with 
the  scaled  versions  of  relative  position  and  velocity 


a 

b 

P 

I 

<1 

e 


^max 

(P 


2yo  +  ^nxo  1  •  .  . 

- =  -Vo  +  4:Xo 

n  TT 

ntjo  -2xo  ^  1  • 

- =  yo - Xo 

n  TT 


\l  + (t)  + (s®”) 

ZgXo  —  Zo'n?{a  —  Xo)  z^Xo  —  4:Tr^Zo(a  —  Xo) 

xl  +  n'^{a- Xo^  il  +  47v‘^{a  -  XoY 

ZoXoTi  +  Zon{a  —  Xo)  t^ZoXo  +  7rZo(a  —  Xo) 

2[xl  +  n^{a- Xo)^]  +  4:TT^{a  -  Xo)^ 


tan  ^ 

n{xo  —  a) 

=  tan  ^ 

27r(5;o  —  a) 

Xo 

Xo 

=  tan  ^ 


(98a) 

(98b) 

(98c) 

(98d) 

(98e) 

(98f) 

(98g) 

(98h) 


and  the  relative  velocities 


•  ~  Z(7l  ~ 

x{t)  =  — x{t)  =  27rpcos(27rt  +  9) 

Th 

•  ~  Stt  ~ 

y(i)  =  — y(t)  =  — 47rpsin(27rt  +  9)  —  Sna 
n 

•  ^  Stt  ~  ~ 

z{i)  =  — z{t)  =  27r/pcos(27rt  +  9)  —  47rgpsin(27rt  +  9) 


or 


•  ^  Stt  ~ 

z{i)  =  — z{t)  =  — 27rZmax  sin(27rt  +  0) 
n 


G,  1  Units  Analysis 

The  variables  x,  y,  and  ^  are  positions  and  defined  in  terms  of  a  length.  Like¬ 
wise,  X,  y,  and  z  are  velocities  and  are  expressed  as  Finally,  n  is  the  angular 

frequency  of  the  orbit  and  is  expressed  in  units  of  These  units  will  not  change 


143 


if  calculated  using  relative  position  and  velocities  as  functions  of  chief  orbit  fractions. 


a  = 


+  ^length 

timFi  timfi  ^  timFi 


length 


time 


time 


b  = 


time 


length  +  ^ 


length 

=  length 


time 


time 


p  = 


\ 


/  length  \ 


(length  —  length)^  +  (  j  =  y  length^  +  length^  =  length 

\  time  / 


I  = 


Q  = 


‘Zl" S-  -  length  (^f  (length  -  length)  ('^Y 
(Zf)"  +  (tZ)"  (length  -  length)" 

length‘^P^*^P^ — h  ]_  _  length) 

®  time  time  time  time  \  o  o  / 


length^ 

time^ 


/  length  \  ^  _l_  length^ 
\  time  /  '  time^ 


length^  _|_  length^ 


time^ 


time^ 


+  (jZ)  (length  -  length)" 


/  length  \  ^  length^ 
\  time  )  '  time^ 

/  length  \ 


=  unitless 


=  unitless 


9  =  tan 


-1 


. (length  —  length)  \  ,  /  \ 

time^ - & — ^ ^  ^  |  ^  radiaus 


length 

time 


■) 


I  length  I 
\  time  / 
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Appendix  H.  Transformations  between  the  Inertial  and  Relative 

Reference  Frames 

In  order  to  use  an  inertially  propagated  truth  model,  we  need  to  be  able  to 
transform  position  and  velocity  vectors  from  the  relative  to  inertial  frames  and  vice 
versa. 


Referencing  Figure  88,  the  position  vector  is 

IJKjR  —  XYZp  _|_  ijkT 


The  velocity  is  calculated  via  the  inertial  derivative  (see  Appendix  A. 2) 


UK 


M  = 


X 


XYZ 


P 


UK 


L 


The  vectors  expressed  in  the  relative  frame  need  to  be  transformed  to  the  inertial 
frame.  If  inertial  position  and  velocity  vectors  of  the  chief  satellite  are  available,  the 
rotation  matrices  can  be  calculated  as 


^R2I 


a: 


Y 


Z 
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where  X,  Y,  and  Z  are  unit  vectors  for  each  direction  expressed  in  the  UK  frame. 
The  superscript  R2I  indicates  the  rotation  matrix  from  the  relative  to  inertial  frame 
and  I2R  indicates  the  rotation  matrix  from  the  inertial  to  relative  frame.  The  radial 
direction  is  simply  the  unit  vector  of  the  chief’s  position 


X  = 


ijKf^ 


If  the  chief  is  in  a  circular  orbit,  then  its  velocity  vector  will  always  be  perpendicular 
to  the  position  vector 

and  the  out-of-plane  direction  completes  the  right-handed  coordinate  system 


iJKf 

Z  =  X  xY  =  X 


Thus,  the  position  vector  can  be  expressed  as 


=  C 


R2I 


and  the  velocity  vector  as 


UK 


M  =  C 


=  C 


UK 


L 


'R2I 


X 

0 

X 

\ 

y 

+ 

0 

X 

y 

z 

n 

z 

/ 

^uk-r 


'R2I 


X  —  ny 
ij  -\-  nx 
z 


+  '■"'X 


(100) 


(101) 
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Equations  (100)  and  (101)  can  be  solved  for  the  relative  vectors  if  we  desire  to  go 
from  the  inertial  to  the  relative  frame. 


C 


R2I 


IJKj 


(102) 


X 

-ny 

(^R2I 

y 

z 

nx 

0 

Noting  that  for  orthonormal  rotation  matrices  [24] 

(^I2R  ^  ^R2I  _  ^I2R  ^  ^  ^I2R  ^  (C''2R^-1 


we  get 


XYZ 


P  = 


=  C^ 


-I2R  /  UK 


M  -  '■""£ 


X 

ny 

XYZp  ^ 

y 

z 

^  ^I2R  (^ijk^  _  ijkyyj  ^ 

—nx 

0 

(103) 


(104) 


An  alternative  formulation  of  the  rotation  matrix  is  derived  henceforth.  By 
inspection,  the  rotation  matrices  are 


X 

cos  i  sin  i  0 

i 

Y 

= 

—  sin  i  cos  £  0 

j 

Z 

1 

O 

O 

k 

CK21 


(105) 
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(106) 


/ 

cost' 

—  sin  £  0 

X 

J 

= 

sint 

cost  0 

Y 

k 

0 

0  1 

y 

Z 

CI2R 


where  £  for  the  circular  orbit  is 


£  =  nto 


and  to  is  the  time  since  passage  over  /.  Equations  (103)  and  (104)  can  be  combined 
into  a  more  compact  form.  Let 


then 


C-I2R 


Xd-Xc 
Yd-Yc 
Zd  —  Zq 


Xc 

4 

1 _ 

= 

= 

Yd 

it; 

Yd 

_ 1 

_Zd_ 

_ 1 

that 

ny 

0  1  0 

X 

—nx 

=  n 

-10  0 

y 

0 

0  0  0 

z 

z 

X 

y 

z 


C12R 


Xd-Xc 
Yd-Yc 
Zd  —  Zc 


+  n 


0 

-1 

0 


1  0 

Xd-Xc 

0  0 

C-I2R 

Yd-Yc 

0  0 

1 

_ 1 
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and  substituting  the  rotation  matrices,  Equation  (106),  we  get 


X 

cos£ 

—  sint' 

0 

0 

0 

0 

-Xc 

y 

sin£ 

cost' 

0 

0 

0 

0 

Yd 

-Yc 

z 

0 

0 

1 

0 

0 

0 

Zd 

-Zc 

X 

nsin£ 

ncost 

0 

cost 

—  sint 

0 

Xd 

-Xc 

y 

—n  cos  £ 

nsm£ 

0 

sint 

cost 

0 

Yd 

-Yc 

z 

0 

0 

0 

0 

0 

1 

Zd 

-Zc_ 

(107) 
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Appendix  I.  Circular  Chief  Relative  Velocities 


Assuming  we  have  found  points  of  interest  on  or  within  the  lobe  between  which 
we  want  the  deputy  to  travel,  a  solution  for  the  initial  and  final  velocities  is  required. 


I.l  Initial  Relative  Velocity 

The  following  is  a  modification  of  the  derivation  presented  by  Mullins  [37]. 
Assume  that  the  deputy  is  at  the  initial  point  at  to  =  0  and  at  the  final  point  at 
tf  =  T.  Let  x{T)  =  Xf,  y{T)  =  y/,  and  z{T)  =  Zf.  Start  with  a  slightly  modified 
form  of  Equation  (82)  from  Appendix  G 


Xf 

Vf 

c 


x(T)  =  —  sin(nT)  - 
n 


2^0 


n 


3Xn 


‘2y 

cos(nT)  H - -  +  Axc 

n 


=  !/(T)=2 


2% 


n 


+  3xo 


2t 

sin  (nT)  H - -  cos  (nT) 

in 


-n 


2^0 


n 


+  dXo 


T  +  Vo 


2Xn 


n 


Zo 


ziT)  =  —  sin(nT)  +  Zq  cos(nT) 
n 


Combining  terms  and  preparing  for  matrix  form 


1  2 

Xf  =  —  sm{nT)xo  H —  [1  —  cos(nr)]  ^o  +  [4  —  3  cos(nr)]  Xg 
2  1 

yf  =  —  [— 1  +  cos(nr)]  Xo-\ —  [4  sm{nT)  —  3nT]  ijo  +  [6  sin(nr)  —  6nT]  Xo  +  yo 
Zf  =  —  sm{nT)zo  +  cos{nT)zo 


yielding  {S  =  sin(nT),  C  =  cos{nT)) 


Xf 

-  [1  -  C'j  0 

n  n  L  J 

Xo 

[4  -  3(7]  Xo 

Vf 

= 

q-l  +  C]  i[4S-3nr]  0 

Vo 

+ 

[bA  -  6nT]  Xo  +  yo 

c. 

0  0  is 

n  _ 

- 1 

O 

(108) 


Note  that  there  is  a  sine  term  in  the  3-by-3  element  of  the  first  matrix.  When  this 
sine  function  goes  to  zero,  the  matrix  is  singular.  Physically  the  deputy  satellite  is 
passing  through  the  Y  axis  (or  a  line  parallel  to  it),  the  intersect  point  for  an  infinite 
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number  of  relative  orbits.  Therefore  the  initial/final  conditions  and  time-of-flight  do 
not  uniquely  define  an  orbit-thus  an  indeterminate  matrix.  This  restriction  should 
not  be  a  problem.  Multiplying  each  side  by  n  and  rearranging  terms: 


nxf  —  n  [4  —  3C]  Xg 

S  2-2C  0 

Xo 

nyf  —  n  [OS'  —  6nT]  Xo  —  nyo 

= 

-2  +  2(7  4,5 -3nT  0 

Vo 

nzf  —  nCzo 

1 

0 

0 

_ 1 

Zo 

Finding  the  inverse  via  co- factors 

=  (4S -3nT)S -0  =  4S^ -3nTS 
Au  =  -l(-2  +  2C)S-0]  =  2S-2SC 
yli3  =  0-0  =  0 

A21  =  -1(2  -  2C)S  -  0]  = -2S  +  2SC 
^22  =  S*S-0  =  S^ 

A23  =  -[0-0]=0 
ylgi  =  0-0  =  0 

A32  =  -[0-0]=0 

^33  =  S(4S-3nT)-(-2  +  2C)(2-2C)  =  4S^-3nTS  +  4-8C  +  4C^ 
=  8-3nTS-8C 

The  determinant  can  be  found  from  the  first  row 

ll^ll  =  +  ^^12^12  +  ^^13^13 

=  ,5(4,52  -  3nTS)  +  (2  -  2C)  (2S  -  2SC)  +  0*0 
=  S{8-3nTS  -8C) 

and  finally  the  inverse  [53] 

^-1  ^  adj  A 
det  A 
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where  the  adjoint  of  A  is  the  transpose  of  the  co- factor  matrix 


A-i  = 


1 


S(8  -  3nTS  -  8C) 


4S^-3nTS  2S-2SC 


-2S  +  2SC 
0 


52 

0 


0 

0 

3nTS  -  8C 


AS'^-SnTS 


-2S+2SC 


S{8-3nTS-8C)  5(8-3nT5-8C) 

2S-2SC 

S{8-3nTS-8C)  S{8-3nTS-8C) 

0  0 


A-'  = 


0 
0 

8-3nTS-8C 
S(8-3nTS-8C)  _ 

—2+2  cos(nT) 

S—3nT  sin(nT) — 8  cos(nT)  8— 3nT  sin(nT) — 8  cos(nT) 
2— 2cos(nT)  sin(nT) 

8—3nT  sin(nT) — 8  cos(nT)  8—3nT  sin(nT) — 8  cos(nT) 

0  0 


4sin(nT)— 3nT 


0 

0 

1 

sin(nT)  _ 


Therefore  the  solution  to  the  initial  relative  velocities  is 


Xo 

Vo 

z, 


4  sm{nT)—3nT 

— 2+2  cos(?T-+') 

0 

8—3nT  sin(nT)— 8  cos(nT) 

8—3nT  sin(nT)— 8  cos(nT) 

2—2  cos(nT) 

sin(nT) 

0 

8—3nT  sin(nT) — 8  cos(nT) 

8—3nT  sin(nT) — 8  cos(nT) 

0 

0 

1 

sin(nT)  _ 

nxf  —  n[4:  —  3cos(nT)]  Xg 
nyf  —  n  [6  sm{nT)  —  6nT]  Xo  —  nyo 
nzf  —  ncos{nT)zo 


We  can  further  break  down  the  second  matrix  into 


n 


— 4  +  3cos(nr)  0 

— 6sin(nT)  +  6nT  —1 
0  0 


0  10  0 
0  0  10 
cos(nT)  0  0  1 


Xo 

Vo 

Zo 

Xf 

Vf 
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Multiplying  this  out,  the  transformation  matrix  becomes 


n 


[4S'-3nr][-4+3C]  +  [-2+2C][-6S'+6nr]  -[-2+2C] 

8-3nTS-8C  8-3nTS-8C 

[2-2C]  [-4+3C]+S'[-6S'+6nr]  -S 

8-3nTS-8C  8-3nTS-8C 

0  0 


0 

0 


c 

s 


AS-3nT 

8-3nTS-8C 

2-2C 

8-3nTS-8C 

0 


-2+2C 

8-3nTS-8C 

S 

8-3nTS-8C 

0 


Simplifying 


n 


-4.S+3nTC 

8-3nTS-8C 

-14+6nT5+14C 


2-2C 


8-3nTS-8C 

-S 


8-3nTS-8C  8-3nTS-8C 

0  0 


0 

0 


c 

s 


4.S-3nT 

8-3nTS-8C 

2-2C 

8-3nTS-8C 

0 


-2+2C  Q 

0 


8-3nTS-8C 

S 

8-3nTS-8C 

0 


Thus  the  transformation  is 


Xo 

-4S^3nTC 

2-2C 

0 

4S-3nT 

-2+2C 

0 

8-3nTS-8C 

8-3nTS-8C 

8-3nTS-8C 

8-3nTS-8C 

Vo 

=  n 

-14+6nrS'+14C 

-S 

0 

2-2C 

S 

0 

8-3nTS-8C 

8-3nTS-8C 

8-3nTS-8C 

8-3nTS-8C 

Zo 

0 

0 

c 

s 

0 

0 

1 

s_ 

Vo 

Vf 

(109) 

We  can  also  express  this  in  terms  of  fractions  of  a  chief  orbit  period  by  substituting 
nT  =  2'kT  (Equation  (2)  in  Section  3.1) 


Xo 

-4iS+67rfC 

2-2C 

0 

AS—GttT 

-2+2C 

0 

8-67rT5-8C 

8-67rT5-8C 

8-6nfS-8C 

8-67rT5-8C 

Vo 

=  n 

-14+127rf'S+14C 

-S 

0 

2-2C 

S 

0 

8-6nfS-8C 

8-6nfS-8C 

8-6nfS-8C 

8-6nfS-8C 

Zo 

0 

0 

c 

s 

0 

0 

1 

s. 

(110) 

where  T  =  T/P  S  =  sin(27rT')  C  =  cos(27rT') 
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If  we  choose  units  of  relative  velocity  that  are  functions  of  time  as  a  fraction 
of  chief  orbit  radius  as  opposed  to  absolute  time,  we  can  use  the  transformation  in 
Equation  (78),  Appendix  F 


Xo 

271 

Xo 

=  27r 

Vo 

n 

Vo 

-4:S+6nTC 

8-6nfS-8C 
-U+l2TTfS+lAC 

S-eTrfS-SC 

0 


2-2C 


8-&'kTS-8C 

-S 

8-6-KfS-8C 

0 


0 

0 


c 

s 


AS—QttT 

3-67rf5-8C 

2-2C 

^-GttTS-SC 

0 


-2+2C 

-6nfS-8C 

S 

-enfs-sc 

0 


(111) 

Finally,  it  is  interesting  to  note  that  the  Po  and  pf  terms  in  the  time-of-flight  matrix 
are  the  negative  of  the  other,  thus  we  can  substitute  a  Ap. 


Xo 

Vo 

Vf 


Xo 

-45+67rfC 

0 

AS-GttT 

0 

-2+2C 

8-67rf5-8C 

8-6ntS-8C 

8-6nfS-8C 

Vo 

=  271 

-14+127rf'S+14C' 

0 

2-2(7 

0 

S 

8-67rT:S-8C' 

8-67rfS-8(7 

8-6-KfS-8C 

0 

c 

s 

0 

1 

s 

0 

where  Ap  =  pf  -po 


Xo 

Xf 

Ay 


(112) 


1.2  Final  Relative  Velocity 

To  calculate  the  velocities  at  the  end  of  an  arc,  we  start  with  the  derivative  of 
Equation  (82)  from  Appendix  G 


Xr. 


Xf 

Vf 


=  x{T)  =  — ncos(nT) 


n 


n 


+  3xo 


-n)  sin(nT) 


/  m\  2^0  /  \  /  m\  ^ 

ncosinT)  H - (—n)  sminT ) - n 

—  +  4:Xo 

n 

n  2 

n 

=  !/(T)  =  2 


ziT)  =  —ncosinT)  +  Zo{—n)  sminT) 
n 
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Preparing  for  matrix  form 


Xf  =  cos{nT)xo  +  2sm{nT)yo  + ‘insm{nT)xo 

yf  =  —2sm{nT)xo+[—3  +  4:Cos{nT)]yo+[—Qn  +  6ncos{nT)]xo 

Zf  =  cos{nT)zo  —  nsm(nT)zo 


yields  (S  =  sm{nT),  C  =  cos(nT)) 


Xf 

Vf 

= 

C  2S  0 
-2S  -3  +  4(7  0 

0  0  (7 

Substituting  Equation  (109),  the  first  term  is 


Xq 

Vo 

+ 

Zo 

3nS  0  0 

-6n  +  6n(7  0  0 

0  0  -nS 


Xq 

Vo 


(7 

27 

0 

-4S'+3nrC 

2-2C 

0 

AS-SnT 

-2+2C 

0 

8-3nTS-8C 

8-3nTS-8C 

8-3nTS-8C 

8-3nTS-8C 

-27 

-3  +  4(7 

0 

n 

-14+6nrS'+14C 

-S 

0 

2-2C 

S 

0 

8-3nTS-8C 

8-3nTS-8C 

8-3nTS-8C 

8-3nTS-8C 

0 

0 

(7 

0 

0 

c 

s 

0 

0 

1 

Xo 

Vo 

Xf 

Vf 


and  the  second  term 


n 


3S  0  0  0  0  0 

-6  +  6(7  0  0  0  0  0 

0  0-7000 


Xo 

Vo 

Xf 

Vf 
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Multiplying  out  the  first  term  and  adding  it  to  the  second 


Ci-4S  +  3nTC)  2S{-U  +  6nTS  +  UC)  3S{8  -  3nTS  -  8C) 

“  “  8-3nTS-8C  ^  8  -  3nTS  -  8C  ^  8  -  3nTS  -  8C 

-AS + 3nT 

8  -  3nTS  -  8C 

^  C{2-2C)  2S{-S)  -2  +  2C 

8  -  3nTS  -  8(7  8  -  3nTS  -  8(7  8  -  3nTS  -  8(7 

Ai3  =  0  +  0  =  0 

C{AS  -  3nT)  2S{2  -  2(7)  AS-  3nTC 

“  8  -  3nTS  -  8(7  ^  8  -  3nTS  -  8(7  ^  “  8  -  3nTS  -  8(7 
C{-2  +  2(7)  25(5)  2-2(7 

“  8  -  3nTS  -  8(7  ^  8  -  3nTS  -  8(7  ^  ^  8  -  3nTS  -  8(7 

Aiq  =  0  +  0  =  0 


-2S{-AS  +  3nTC)  (-3  +  4(7)(-14  +  6nT5  +  14(7)  (-6  +  6(7)(8  -  3nr5  -  8(7) 

“  8  -  3nTS  -  8(7  ^  8  -  3nTS  -  8(7  ^  8  -  3nr5  -  8(7 

2-2(7 

8  -  3nr5  -  8(7 

^  -25(2 -2C)  (-3  +  4C)(-5)  -5 

8  -  3nTS  -  8(7  8  -  3nr5  -  8(7  8  -  3nTS  -  8(7 

^23  =  0  +  0  =  0 

-25(45 -3nr)  (-3  +  4(7) (2  -  2(7)  -14  +  6nT5  +  14(7 

“  8  -  3nTS  -  8(7  ^  8  -  3nTS  -  8(7  ^  “  8  -  3nr5  -  8(7 

^  -25(-2  +  2C)  (-3  +  4C)(5)  5 

8  -  3nTS  -  8(7  8  -  3nT5  -  8(7  8  -  3nTS  -  8(7 

^26  ~  0  +  0  =  0 


A^i  —  0  +  0  —  0 

As2  —  0  +  0  =  0 

^33  =  C^  +  (-S)f  =  -i 

Asa  —  0  +  0  =  0 
^35  =  0  +  0  =  0 

-4*  =  ei  +  o=2 
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Thus  the  relative  velocity  at  the  end  of  a  trajectory  can  be  expressed  as 


Xf 

-4S+3nT 

-2+2C 

0 

4S-3nTC 

2-2C 

0 

8-3nTS-8C 

S-SnTS-SC 

8-3nTS-8C 

8-3nTS-8C 

Vf 

=  n 

2-2C 

-S 

0 

-14+6nrS'+14C 

S 

0 

8-3nTS-8C 

8-3nTS-8C 

8-3nTS-8C 

8-3nTS-8C 

0 

0 

1 

s 

0 

0 

c 

s_ 

(113) 

We  can  also  express  this  in  terms  of  fractions  of  a  chief  orbit  period  by  substituting 
nT  =  2'kT  (Equation  (2)  in  Section  3.1) 


Xo 


Xf 

-45+67rT 

-2+2C 

0 

4S-67rfC 

2-2C 

0 

8-67rTS'-8C 

8-6nfS-8C 

8-6nfS-8C 

8-67rT5-8C 

Vf 

=  n 

to 

1 

to 

-S 

0 

-14+127rf5+14C 

S 

0 

8-67rTS'-8C 

8-6nfS-8C 

8-67rT5-8C 

8-6nfS-8C 

0 

0 

1 

s 

0 

0 

c 

s. 

where  T  =  T/P  S  =  sin(27rT')  C  =  cos(27rT') 


( 


14) 


If  we  choose  units  of  relative  velocity  that  are  functions  of  time  as  a  fraction 
of  chief  orbit  radius  as  opposed  to  absolute  time,  we  can  use  the  transformation  in 
Equation  (78),  Appendix  F 


Xf 

-45+67rf 

-2+2C 

0 

4S-67rfC 

2-2C 

0 

27r 

8-67rf5-8C 

8-6nfS-8C 

8-6nfS-8C 

8-6nfS-8C 

if 

Vf 

=  27r 

2-2(5 

-S 

0 

-14+127rT:5+14(5 

S 

0 

n 

S-eirfS-SC 

S-eirfS-SC 

8-67rfS-8(5 

8-6TTfS-8C 

0 

0 

1 

s 

0 

0 

c 

s. 

(115) 
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Similar  to  the  initial  relative  velocity,  we  can  express  this  in  terms  of  Ay 


Xf 

-4*§+67rf 

0 

AS-^ttTC 

0 

2-2C 

8-67rfS-8C 

8-67rfS-8C 

8-67rfS-8C 

h 

=  27r 

2-2C 

0 

-U+UttTS+UC 

0 

S 

8-67rT5-8C 

8-6nfS-8C 

8-6nfS-8C 

0 

1 

s 

0 

c 

s 

0 

(116) 


where  Ay  =  yf  -po 


158 


Appendix  J.  Tmax  of  the  Z  Motion  for  a  Cireular  Chief 

The  Z  motion  is  decoupled  from  the  XY  motion  and  can  thus  be  treated  sepa¬ 
rately.  The  linearized  equation  of  motion  is 

2  -|-  n^z  =  0 

which  is  a  harmonic  oscillator.  The  solution  to  this  differential  equation  is  (Equation 
(94),  Appendix  G) 


z{t)  = 


nt  —  tan 


-1 


nzr 


The  maximum  and  minimum  .2  values  are  defined  as  follows 


^max  =  furthermost  allowable  z  coordinate  from  the  chief 
^min  =  closest  allowable  z  coordinate  to  the  chief 

Pulling  the  amplitude  information  from  the  above  equation  and  setting  it  equal  to 


=  n^/. 


^max  - 


or 


nzo 


(117) 


Assume  a  cylinder  is  defined  as  shown  in  Figure  89  where  h  is  half  the  height  of  the 
cylinder  (h  >  0)  and  (3  is  the  angle  between  the  Z  axis  and  the  vector  from  the  chief 
to  the  center  of  the  lobe.  We  wish  to  find  the  maximum  time-of-fiight  (Tmax)  between 
a  starting  position  Zo  and  final  position  z/.  Assume  that  the  deputy  starts  and  stops 
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(b)  f  <  /?  <  TT 


Figure  89:  Cylindrical  Lobe 


at  the  same  position,  then 


Zf  =  Zo 


COS 


TiTyY^pi,-^  tan 


-1 


Solving  for 


cos 


-1 


.•JW+ 


COS 


-1 


^  +1 
nzo  ' 


tan 


-1 


nzn 


77/^7Tnov  tan 


=  tiT 

-  I  L±  YY 


nZo 


cos 


-1 


^  +1 

nzo  ' 


+  COS 


-1 


1 


+1 

Y  nzo  ' 


=  nT„ 
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Collecting  terms  and  substituting  Equation  (117) 


Tmax  COS 

n 


-1 


^  +1 

nzc 


2 

—cos 
n 


-1 


\ 


2^0 


-1+1 


=  —cos 
n 


-1 


Assume  that  the  deputy  enters  the  lobe  at  =  ^min  and  then,  by  Figure  89, 


=  ^min  =  J  COS  P  -  eh 
^max  T  COS  /3  “h  sH 


(118) 

(119) 


where 


Thus 


where 


£  = 


1  if  0  <  ^ 

-1  if  —  <  /3  <  jr 
2  ~ 


2  f -f  cos /3  -  eh\  2  _i  f  z-^in 

Jmax  =  -COS  - — - -  =  -COS  - 

n  V7COS+  +  £hy  n  V^max 


0  < 


7  cos  13  —  eh 


<  1 


(120) 


(121) 


7  cos  (3  -\-  eh 

The  lower  bound  occurs  when  the  cylinder  is  tangent  to  the  orbital  plane  while  the 
upper  bound  occurs  when  the  cylinder  has  zero  height.  The  range  of  the  inverse 
cosine  on  the  domain  [0, 1]  is  [O,  |] ,  which  means  the  domain  of  time-of-flight  is  [0, 1] 
(reference  Figure  90). 

These  can  also  be  written  as  functions  of  chief  period  in  which  the  linearized 
equation  of  motion  is 

z  +  dyr^z  =  0 
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FigUrG  90.  Tmax  ^min/^max 


yielding  a  closed  form  solution  of 


z{i)  =  —  sin(27rf)  +  5^  cos(27rt) 
2ti 


Applying  the  harmonic  addition  theorem  (Appendix  A. 3) 


z{i)  = 


27rt  —  tan 


-1 


2'KZr 


thus 


and  the  maximum  time-of-flight  is 


Tmax  1  _i  /7COS/?  - 
=  — —  =  —COS  - - - - 

P  7T  \7c0s  p  +  sh  J 


—  COS 
TT 


(122) 


(123) 
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Appendix  K.  The  Derivation  of  AV 

This  research  is  founded  on  the  assumption  that  the  deputy  satellite  can  perform 
impulsive  burns  to  keep  itself  inside  the  defined  lobe.  Calculating  the  magnitude  of 
those  impulsive  burns  is  therefore  a  key  component  of  the  cost  function  used  in  the 
optimization  process.  The  following  derivation  finds  the  magnitude  of  the  impulsive 
burn  (AC). 

K,  1  The  Derivation  of  AC^ 


Figure  91:  AC  Derivation  Setup 

Avi  is  defined  as  the  difference  between  the  velocity  before  and  after  an  impulsive 
thrust,  that  is,  the  instantaneous  change  in  velocity  required  at  the  burn  point  [62] 
and  is  shown  vectorially  in  Figure  91. 

Avi  =  ^1  -V- 

The  Z  component  completely  decouples  from  the  XY  thus  each  will  be  addressed 
separately.  In  the  orbit  plane  of  the  chief 

vt  =  xtX  +  yfY 
C  =  + 
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Let  the  magnitude  of  the  Av  vector  be 


=  ||Aui||2  = 

The  mathematics  of  the  derivatives  simplify  somewhat  if  we  remove  the  square  root 
and  deal  with  Al/^.  This  will  not  adversely  affect  our  ability  to  use  AV  in  a  cost 
function.  Since  Al^  is  always  a  positive  value,  minimizing  the  square  is  equivalent  to 
minimizing  the  square  root. 


At-2 


(4 


-^2 


)  +  (%  -  Vi  ) 


x'l  X- 

vt  -  Vi 


(124) 


where  the  velocity  prior  to  the  burn  (Equation  (113),  Appendix  I)  is 
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where  T  =  S  =  sm(nTi-i^i),  C  =  cos(nTj_i^j)  and  the  velocity  after  the 

burn  (Equation  (109),  Appendix  I) 


4 

0 

0 

-4S++3nT+C+ 

2-2C+ 

4S+-3nT+ 

-2+2C+ 

=  n 

8-3nT+S+-8C+ 

8-3nT+S+-8C+ 

8-3nT+S+-8C+ 

8-3nT+5+-8C+ 

Vi 

0 

0 

-l4+%nT+ S+ +14C+ 

-S+ 

2-2C+ 

S+ 

8-3nT+S+-8C+ 

8-3nT+S+-8C+ 

8-3nT+S+-8C+ 

8-3nT+S+-8C+ 

where  T+  =  =  sm{nTi^i+i),  =  cos(nTj,j+i).  Let 


Xi-i 

Vi-i 

Xi 

Vi 

^i+l 

Vi+I 


(125) 
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Vi  -  Vi 


=  n 


^i—l  Vi—l  Ui  ^i+1  Vi+l 


R 


Then 


R  = 


45--3nT- 

8-3nT-S--8C- 

2-2C- 

8-3nT-  S-  -8C- 

-4S++3nT+C+  _  4S--3nT-C 

8-3nT+S+-8C+ 


2-2C+ 


8-3nT-S--8C- 

_ 2-2C- 

8-3nT+S'+-8C+  8-3nT-S--8C- 

4S+-3nT+ 

8-3nT+S+-8C+ 

-2+2C+ 

8-3nT+S+-8C+ 


-2+2C- 


i-3nT- S-  -8C- 

s- _ 

i-3nT-  S-  -8C- 


-l4+&nT+S++l4C+ 

8-3nT+S+-8C+ 

-S+ 


-l4+6nT- S- +14C- 

8-3nT-S--8C- 


8-3nT+S+-8C+  8-3nT- S- -8C- 

2-2C+ 

8-3nT+S+-8C+ 

_ ^ _ 

8-3nT+5+-8C+ 


(126) 
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We  can  also  express  this  matrix  in  terms  of  fractions  of  a  chief  orbit  period  by  sub¬ 
stituting  nT  =  2tiT  (Equation  (2)  in  Section  3.1) 


R  = 


where 


4S--67rf- 

8-67rT-S--8C'- 

2-2C- 

8-6TTf-S--8C- 

-45++67rT+C+  _  4:S--6nf-C- 

8-6nf+S+-8C+  8-6nf-S--8C- 

2-2(7+ _ 2-2C- 

8-6nf+S+-8C+  8-6nf-S--8C- 

4S+-6TTf+ 

8-67rT+S+-8(7+ 

-2+2(7+ 

8-67rf+S+-8(7+ 


-2+2(7“ 

8-67rT-S--8(7- 

s- _ 

8-67rf-S--8(7- 

-14+1 2wf+  5+ + 14(7+  _  -14+127rf'-5-+14(7~ 
8-67rf+5+-8(7+  8-6nf-S--8C- 

-S+ _ SZ _ 

8-67rT+S'+-8(7+  8-6nf-S--8C- 

2-2(7+ 

8-67rf'+S+-8(7+ 

S+ 

8-67rr+S+-8(7+ 

(127) 


T- 

s- 

C- 


Finally  AV)^  is 


sin(27rf'j_i,j) 

cos(27rf)_i,j) 


T’+  = 

=  sin(27rTj,j+i) 
(7+  =  cos(27rTi,i+i) 


AV^  =  r? 


Vi—l  Vi  ^i+1  Ui+l 


RR' 


Xi-l 

Vi-i 

Xi 

Vi 

^i+l 

Vi+I 


(128) 


where  we  can  substitute  R  for  R  as  desired.  AV  can  also  be  expressed  in  terms  of 
velocity  as  a  function  of  fractions  of  a  chief  orbit  period  (see  Appendix  F).  Starting 
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with  Equation  (124) 


Al/.^  = 


n  ^  ,  n  x._  n  n 

— xj - X-  — yj - 
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This  equation  reduces  to  Equation  (128)  if  we  substitute  Equations  (111)  and  (115) 
but  is  useful  if  the  pre-  and  post-scaled  relative  velocities  are  known. 


K.2  Specific  AV 

It  is  useful  to  find  a  form  of  AV  that  is  not  dependent  upon  the  chief’s  semi¬ 
major  axis  (or  orbital  radius  for  circular  orbits).  Therefore,  define  specific  delta  V 
(AE)  as 


AE"  = 


AE" 


Xi—l  Vi—l  Xi  Pi  Xj-|-i  Vi+l 


RE! 


xJ  -  x^  yfi  -  ?/j 


vT  -  Vi 


Xi-l 

Vi-I 

Xi 

Vi 

Xi-\-l 

Vi+i 


(129) 


and  the  conversions  between  different  types  of  AE  are 


AE:  = 


AE  _  AE 

n  2tt 


(130) 
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K.3  The  Derivation  of  AV  for  Z  Motion 

Since  the  Z  component  of  the  deputy’s  motion  decouples  from  XY  plane,  it  can 
be  treated  separately.  The  instantaneous  change  in  velocity  at  the  burn  points  is  [62] 

AV,  =  i+  -  i- 


The  relative  velocity  prior  to  the  burn  (zf)  is  the  final  relative  velocity  from  the  prior 
leg,  and  the  relative  velocity  after  the  burn  (zf )  is  the  initial  relative  velocity  of  the 
next  leg.  These  are  defined  by  Equations  (109)  and  (113) 


Zn  =  -nZn 


COS  (nTj) 

sin(nrj) 

1 


Zf  =  —nZo— 


sin  (nTj) 


+  nzf 
+  nzf 


1 

sin  (nTj) 
cos(nTj) 
sin(nTj) 


Thus 


A14,  = 


cos  (nTd  1 

-nZn— — .  +  nzf  — 


sm{nTi)  "sin(nr,)_ 


-nzr. 


1  cos(nTj) 

sin(nTj)  sin(nTj) 


=  nZn 


1  —  cos(nTi) 

1  'y  n 

1  —  cos(nTj) 

W 

1 _ 

-h  nzf 

sm{nTi) 

=  n{zo  +  Zf)  tan 

or  in  terms  of  chief  orbit  fractions  by  substituting  Equation  (2),  Section  3.1 


AVj  =  n{zo  +  Zf)  tan 


Note  that  AE  is  linear  in  Zg  and  Zf;  thus  choosing 


'2'o  ^min 
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will  result  in  the  minimum  AV 


AVi  =  2n2;mm  tan  (131) 

This  substitution  also  allows  the  following  reduction  in  the  relative  Z  direction  veloc¬ 
ities 


•  _  l-cos(nTi) _ 

Zq  ^^min  •  /  rri  \  ^^min  tail  I  TiT^ 
sm[nTi)  \2 

^  l-cos(nri)_ 

sin(nT^)  \2 


or  in  terms  of  chief  orbit  period 


Z(^  27r^T-nin  tan  (  tvT^ 


Zf  =  —  27rimin  tan  (  ttT; 


Assume  the  total  hover  time  (Tt)  is  given,  the  period  of  the  Z  motion  is  constant 
(Pz),  and  that  the  entry  conditions  for  the  Z  motion  is  open,  then  the  total  number 
of  burns  required  to  keep  the  the  deputy  inside  the  lobe  is 


^  Burns 


and  the  AV  required 


AVz  —  2?r 


Zr^^in  tan  (  ttP; 


or,  in  terms  of  specific  AV,  (see  Appendix  K.2) 


(132) 


AVz 

AVz  =  — ^  =  2 


n 


Tt 

Pz 


Zry,\y,  tan  (  TlP, 


(133) 


where  [  J  represents  the  floor  function.  Equation  (132)  converges  to  the  continuous 
AV  solution  (Equation  (137)  in  Appendix  L)  as  Pz  goes  to  zero  and  is  proven  below. 
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Let  Pz  0.  Using  the  small  angle  approximation, 


tan  (  ttPz  1  = 


sin  (  ttP, 


^  TvPz 


cos  (  ttP, 


=  ttPz 


Further,  let  Pz  be  an  whole  fraction  of  TV  where  kz  is  the  number  of  Z  legs 


^Tt  =  P, 


Then 


AUz  = 


Ik. 


■Tt 


*  2nz^in  (  tt^Tt 
kz 


=  2n^min7r  [kz\  t^Tt 

Kjz 


Since  kz  is  an  integer 


thus 


or  in  terms  of  specific  AU 


[kz\  =  kz 


AVz  =  2nZn,in7rTT 


AW 

AW  = - = 

n 


which  is  the  same  as  the  continuous  burn  solution  (Equation  (137)  in  Appendix  L). 
The  constraint  on  choosing  kz  is 


Z—Tt  —  Pz  <  Umax 
kz 


Solving  for  kz 


kz  > 


Tj 

Th. 
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since  is  an  integer 


(134) 

where  [  ]  represents  the  ceiling  function.  The  total  AV  expended  performing  these 
optimal  burns  is  therefore 


k.  = 


^Vz  =  {kz  -  1)  *  2n2(mintan  (  j  =  {k^  -  1)  tan  |  ^ 


in  specihc  AV 


AVz  =  25min  {kz  -  1)  tan 


ttTt 


or  expressed  in  terms  of  the  initial  or  final  velocities  in  the  Z  direction 


AVz  =  2{kz  -  l)zo  =  -2{kz  -  l)zf 


and  in  terms  of  velocities  that  are  functions  of  the  chief’s  period  and  specific  AV 


•  77  •  77 

A14  =  2{kz  -  l)5o^  =  -2{kz  -  1)%  — 


Thus 


AVz  =  -{kz  -  l)zo  =  --{kz  -  l)zf  (135) 

TT  TT 

Figure  92  represents  an  example  plot  of  Equation  (133).  Note  the  sawtooth  feature 
that  is  an  artifact  of  the  floor  function  and  that  the  discrete  solution  converges  to  the 
continuous  solution  as  ^  0. 


K.4  AV  Expressed  in  the  Inertial  Frame 

AV  can  also  be  expressed  in  inertial  frame  via 


AE  = 
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Figure  92:  Notional  AV  vs  Z  Period 


Substituting  Equation  (101) 
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Appendix  L.  The  Continuous- Thrust  Solution 

As  a  benchmark  with  which  to  compare  our  optimal  trajectories,  the  continuous- 
thrust  solution  is  derived.  We  can  quickly  find  the  closed-form  solution  for  a  continuous- 
thrust  controller  attempting  to  keep  the  deputy  at  a  specified  point  in  the  relative 
frame.  Again,  we  expect  the  impulsive  burn  to  be  fuel-optimal  compared  to  the 
continuous-thrust  solution  derived  below.  We  want  the  velocities  at  this  hover  point 
to  be  zero 

5  =  0  y  =  0  z  =  0 

which  reduces  Equation  (81),  Appendix  F  to 

X  =  47r(0)  -|-  =  127r^5; 

y  =  — 47r(0)  =  0 

5  =  — 47r^5 

Since  the  y  equation  is  now  zero,  we  need  only  worry  about  the  accelerations  in  the 
X  and  Z  directions.  Integrating  from  zero  to  the  total  time-of- flight  (Tt)  yields  the 
AV  as  a  function  of  chief  orbit  period  (AW)  required  to  keep  the  deputy  hovering  at 
a  given  Xo,  ijo,  Zo'. 

nTj'  rTj'  rTj'  rTj' 

AVc  =  /  \x\di -\-  /  \z\di  =  /  127r^|xo|dt -I-  /  ATS‘\zo\di 

Jo  Jo  Jo  Jo 

=  11“^  {12\x o\  +  4:\zo\)  f  di  =  {12\xo\-\-4:\zo\)fT 

Jo 

where  the  subscript  C  represents  “continuous” .  To  use  this  formulation,  we  must  be 
careful  not  to  pick  lobe  shapes  that  recouple  the  in-plane  and  out-of-plane  compo¬ 
nents.  Note  that  no  claim  is  made  that  this  is  the  optimal  continuous-thrust  solution; 
it  will  merely  serve  as  a  benchmark  controller  that  is  easily  synthesized.  The  location 
of  Xo  and  Zo  is  arbitrary  but  the  smallest  continuous  AV  is  attained  when  Xo  and 
Zo  are  at  their  minimum,  thus  Xmin  and  5min  represent  the  point  on  the  lobe  that  is 
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closest  to  the  Y  axis  and  XY  plane  respectively: 


AVc  = 


(136) 


If  the  lobe  happens  to  intersect  the  Y  axis,  then  the  optimal  solntion  (for  the  linear 
CW  case)  is  to  stay  on  that  Y  axis  which  is  in  fact  the  loci  of  eqnilibrinm  points  and 
requires  no  expenditure  of  fuel.  We  can  also  express  the  continuous  AV  in  terms  of 
specific  AV  by  applying  the  conversions  in  Equation  (130),  Appendix  K 


AW  = 


1 

2tt 


(137) 


The  above  equation  works  well  if  we  assume  the  deputy  starts  at  the  minimum  x  and 
values;  however,  in  order  to  use  this  as  a  fair  comparison,  the  AV  needed  to  get 
into  that  position  should  also  be  included.  Assume  that  the  continuous-thrust  solution 
starts  at  the  same  entry  position  and  velocity  as  the  discrete-thrust  solution.  The  AV 
required  to  get  to  the  Xmin  position  is  AV^  is  (from  Equation  (124)  in  Appendix  K) 


)'  +  -  th  )']  =  hi 


^i  -  yJ-  Vi 


—  Xi 

y{ 

-yi  _ 

Substituting  Equation  (110)  in  Appendix  I 
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• 

yi 
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-14+127rfS+14C' 
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-6nfS-8C 
S 


S-e-n-TSSC  S-e-n-TSSC, 


-1 

Xi 

r  •  1 

X2 

Vi 

Ay 

where 


Ay  =  y2-  yi 
S  =  sin(27rTi^2) 
C  =  cos(27rTi^2) 
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Once  at  the  Xmin  position,  the  deputy  must  make  a  burn  to  cancel  its  relative  velocity 
that  is  equal  in  magnitude  but  opposite  in  direction  to  the  final  velocity  of  the  previous 
leg 


where 


^2 

-ASXGttT 

AS-GttTC 

2-2(7 

_ 

8-6iTfS-8C 

S-67rfS-SC 

S-67rfS-8C 
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Appendix  M.  The  Persistent  Hover  Orbit 

The  persistent  hover  orbit  is  a  special  type  of  repeating  hover  orbit  that  takes 
advantage  of  the  equilibrium  condition  of  the  linear  Clohessy-Wiltshire  equations  of 
motion.  As  presented  in  Section  3.2.2,  \etx  =  y  =  z  =  x  =  y  =  z  =  0.  Then  the 
Clohessy-Wiltshire  equations  reduce  to 

—  127r^x  =  0 
471^5  =  0 

Thus  X  =  0  and  z  =  0  without  any  restriction  on  y,  meaning  the  entire  Y  axis  is 
an  equilibrium  solution.  If  the  user-defined  lobe  is  tangent  to  or  contains  within  it 
any  part  of  the  Y  axis,  then  this  becomes  an  ideal  place  for  the  deputy  satellite  to 
hover.  In  fact,  finding  a  minimum- fuel  solution  to  arrive  and  then  remain  on  either 
the  Y  axis  itself  (single  point  hovering  orbit)  or  any  2x1  ellipse  dictated  by  Equation 
(87),  Appendix  G,  that  is  contained  within  the  lobe.  In  the  absence  of  linearization 
error  and  perturbations,  the  deputy  could  stay  in  this  relative  orbit  indefinitely.  Even 
with  such  nonlinearities  present,  it  would  take  relatively  little  fuel  to  keep  the  deputy 
in  this  persistent  orbit.  Calculating  those  correction  burns  is  not  addressed  in  this 
research. 

When  investigating  these  unique  persistently  hovering  orbits,  two  questions 
must  be  answered:  1)  How  do  we  define  such  a  persistent  orbit,  and  2)  How  do 
we  know  a  candidate  persistent  orbit  fits  within  the  lobe?  The  first  question  is  rela¬ 
tively  easy.  Since  the  eccentricity  of  the  hover  orbit  is  required  to  be  ^  and  it  must 
be  centered  along  the  Y  axis,  only  the  center’s  y  coordinate  (6)  and  relative  orbit  size 
(p)  are  needed  (reference  Eigure  93).  The  second  question  is  answered  via  a  numerical 
solution.  For  a  given  persistent  hover  orbit,  each  point  is  checked  (to  whatever  reso¬ 
lution  the  user  specifies)  for  breaches  of  the  lobe  boundary.  This  is  accomplished  by 
comparing  the  hover  orbit  radius  (rp)  to  the  lobe’s  radius  measured  from  the  hover 
orbit  center  (rp).  If  at  any  point,  rp  >  rp,  then  the  hover  orbit  is  not  acceptable. 
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XL 


Figure  93:  Persistent  Hover  Orbit  Setup 


First,  a  relationship  between  ijj  and  ip'  is  needed.  By  inspection 

xl  +  rcosip 


where  XL,yL  are  the  coordinates  of  the  lobe  center  and  r  is  the  radius  of  the  lobe 
measured  from  its  center.  With  ip'  in  hand,  the  other  two  radii  are  found  with  simple 
trigonometry.  Finding  rp  is  accomplished  via  Equation  (32),  Appendix  A, 


rp 


(2P)(P) 

\/ cos^  Ip'  +  sin^  Ip' 


2p 

a/  4  cos^  Ip'  +  sin^  Ip' 


and  rp  by  the  Pythagorean  theorem 


rp  =  a/  {xp  +  rcosipy  +  {yp  —  yp  +  rsinipy 


where  r  is  defined  by  Equation  (41),  Appendix  B: 


r  = 


-y/ cos^  {ip  —  rj)  +  sin^  {ip  —  rf) 
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Thus  a  numerical  check  can  be  done  by  sweeping  ijj  between  0  and  27r  to  the  desired 
resolution  and  checking  for  rp  >  tl.  Note  that  there  is  a  singularity  when  ip'  passes 
through  I  and  at  which  points 


xl  =  —r  cos  tjj 


The  lobe  radius  can  be  calculated 


tl  =  Vl  —  b  +  r  sin 

but  the  persistent  orbit  radius  is  dependent  on  ■0'  which  is  ill  dehned.  Fortunately, 
these  two  points  are  co- linear  with  the  semi- major  axis 


rp  =  2p 


M.l  Maximum  p  and  b 

In  order  to  search  numerically  for  candidate  persistent  hover  orbits,  we  need  to 
bound  both  p  and  b.  Referencing  Figure  93,  lobes  potentially  intersect  the  Y  axis  at 
two  points 


2/Upper  =  7xy  sin  a  +  ri'ipu)  sin  xpu 
VLowev  =  IxY  sin  “  +  r(^L)  sin  XpL 

The  largest  possible  p  value  is 

7  2/Upper  1/Lower 

r(^[/)  smilJu  -r{i)L)  sin^i 


(138) 
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and  the  range  of  b  by  inspection  is 


femax  =  yUpper  =  7xy  «  +  r{'tpu)  slu  i^u 
bmin  =  yLower  =  7xy  Sin  «  +  r(?/^L)  sin  V'l 
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Appendix  N.  The  Teardrop  Hover  Orbit 


This  appendix  develops  a  proof  that  the  AV  required  to  stay  on  a  teardrop 
orbit  is  equivalent  to  hovering  continuously  at  the  time-averaged  x  coordinate  of 
the  teardrop.  The  results  are  equivalent  to  a  similar  derivation  in  [28].  A  notional 
teardrop  orbit  and  centroid  location  is  shown  in  Figure  94. 
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Figure  94:  Notional  Teardrop  Orbit 


Let  the  subscripts  o  and  /  represent  initial  and  final  conditions  respectively. 
Noting  the  A  scaling  to  convert  to  specific  AV  (Appendix  K),  the  AV  required  for 
each  cycle  of  a  teardrop  orbit  is 

-  Xff  +  (Po  -  Pff  (139) 

Neglecting  motion  in  the  Z  direction,  the  relative  velocities  are  calculated  via  Equa¬ 
tions  (112)  and  (116)  from  Appendix  I. 
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-U^UttTtdS^UC 

2-2(7 

s 

S-GttTtdS-SC 

S-GttTtdS-SC 

S-QttTtdS-SC  _ 

Xf 

Ay 
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if 

-AS^GttTtd 

A^S-Qt^TtbC 

2-2(7 

=  27T 

S-e-KTroS-SC 

^-QttTtdS-^C 

S-GttTtdS-SC 

P. 

2-2(7 

-14+127rfTi:»5+14C 

s 

.8-67rfTDS-8C 

8-67rfTD  5-8(7 

8—67tTtbS—8C_ 

Xo 

Xf 

Ay 


where  Ttd  is  the  time-of-flight  of  a  cycle  (i.e.  period  of  the  teardrop),  S  =  sin(27rTT_D), 
and  C  =  cos{27:Ttd)-  Since  the  initial  and  hnal  coordinates  occur  at  the  apex  of  the 
teardrop, 

Xo  =  Xf 

=  Vo  -  Vf  =  0 

Letting  D  =  8  —  GttTtdS  —  8C, 


Xo 


Vo 


=  2% 


-4S+6ttTtdC 

S-GwTtdS-SC 
— 14-|-127rTVn>S-|-14(7 

.  8-67rfTDS-8C 


4S-6ttTtd 

i-ewfTDS-8C 
2-2(7 

i-6nfTDS-8C2  l'"°2 


if 

-AS+GttTtd 

A^S-Qt^TtbC 

=  27T 

8-&%fTDS-8C 

8-6vfTDS-8C 

yf. 

2-2(7 

-14+127rfTD5+14(7 

_8-67rTTD  5-8(7 

8-6ntTDS-8C 

therefore 

Xo  —  Xf 
Vo-Vf 


27rXo 

2TrXo 

~w 


Xo 

2ttXo 

— 67rTj'£)(l  —  (y) 

- 

Xo 

D 

-12  +  127rrrD«S  +  12(7 
(140) 

-1 

Xo 

27rXo 

Xo 

D 

-12  +  127rfTD^  +  12(7 

(141) 


-67rTh£)(l  —  (5)  —  67rT7’£)(l  —  C*) 


-2A7^^fTDXo{l  -  C) 


D 


-U  +  UttTtdS  +  12C]  -  {-12  +  127vTtdS  +  12C 


=  0 
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Substituting  these  into  Equation  (139)  yields 


f -2ATT^fTDXo{^  -  C) 

\  ^ 


+  (of 


127rrr£)Xo(l  -  C) 

8  —  QttTj'dS  —  SC 


We  can  set  this  equal  to  the  continuous  hover  solution  (Equation  (137),  Appendix  L) 
and  solve  for  the  equivalent  x  coordinate  (xeq)-  This  is  the  x  position  at  which 
thrusting  continuously  would  yield  the  same  AE  as  thursting  at  the  apex  of  the 
teardrop. 


127rTTD5;o(l  —  C) 
8  —  QttTxdS  —  SC 


eqTrTrD 


Solving  for  Xeq, 


2(1  -  C)xo 

8  —  QtvTj'dS  —  SC 


(142) 


This  equivalent  x  coordinate  has  a  physical  meaning.  It  is  the  time-averaged  X 
position  of  the  deputy;  a  proof  is  presented  below.  The  time-averaged  X  position  can 
be  calculated  via 


cTtD 


xtd 


x{i)di 


Ttd  Jo 

where  the  x  position  is  given  by  Equation  (85a)  in  Appendix  G: 


(143) 


x{t) 

x{t) 


ixo  —  a)  cos(nt)  -|-  —  sm(nt)  +  a 
n 


Xn 


{xo  —  a)  cos(27rt)  H - sin(27rt)  -|-  a 

2tt 


and  a  is  Equation  (98a), 


therefore  the  integral  is. 


a  =  —yo  +  4x 

TT 


o 


a)  cos(27rt)dt  -|- 


Jo  27r 


sin(27rt)dt 


adi 
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The  first  term  reduces  to: 


pTtd 

/  {xo  —  a)  co^{27ii)di  = 

Jo 


27T 


1  Ttd 


sin(27rt) 


J  0 


27r 


sin(27rTTD) 


The  second  term  reduces  to: 


rTTD  ^ 

'T/n 


/  —  sm(27ri)di  =  ^ 

/o  27t  ^  ^  47r2 


"I  Ttd 


cos(27rt) 


Xr. 


J  0 


cos(27rTh£))  —  1 


and  the  third  term  reduces  to: 


rTTD 


adt  =  a 


r  1  Ttd 

i 

L  J  0 


—  aTxD 


Putting  it  all  together  and  substituting  for  a, 


rTTD 


x{i)di  = 


~Jo  -  l^io)  S+^^(l-C)+(li„  +  4i,]Trn 


—  (  4Ttd  ~  )  To  H - ^ — —X. 


271 


47|-2 


° ' '  I 


—  3S  _  1  —  C  ji  27rTj^D  —  S  ji 

-Xn  H - ^ — ^Tq  H - ;;; — 7, - Vo 


27T 


471-2 


27r^ 


Substituting  Equation  (140)  for  the  initial  relative  velocities, 


SttTtd-SS 

2tt 


^  (^)  67rT™(l  -C)  +  (^)  -12  +  127rT™^  +  12C 


1 

27tD 


2ttD 


Xn 


SttTj^d  —  35^  j  1?  —  Q7rTj^j7)(^l  —  C)‘^  +  (  27tTxd  —  ^ )  (  — 2A7tTtdS  +  240 
^  IQttTj^d  —  XQitTj'dC  —  —  Q7tTj^]j(^1  —  O)^ 

167rfTD(l  -  O)  -  GTifroiS^  +  (1  -  Cf) 

167tTj^j7}(^1  —  O)  —  127rT7^7^(l  —  O) 

4:7tTj^]j(^1  —  O) 


Xn 


1 

27tD 


Xo 


1 

27ri:> 


Xn 


1 

27rn) 


Xn 
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thus 


pTtd 

Jo 


x{i)di  = 


2TrD(l  ~  C)xo 


—  QttTxdS  —  SC 

The  time-averaged  x  coordinate  is  therefore, 


XtD 


2(1  -  C)xo 

8  —  QtvTj'dS  —  SC 


(144) 


Comparing  this  result  to  Equation  (142),  we  find  that  Xeq  is  in  fact  the  time-averaged 
X  coordinate  and  confirming  that  thursting  continuously  at  xtd  requires  the  same 
amount  of  fuel  as  discretely  thursting  at  the  teardrop  apex  to  remain  on  the  teardrop 
orbit.  The  specific  AV  per  cycle  is 


127rTTD(l  —  C)xo 
8  -  GttTtdS  -  SC 


(145) 
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Appendix  0.  The  Truth  Model 


0.1  Propagation  of  the  Truth  Model 

The  truth  model  is  based  on  the  nonlinear  inertial  orbit  equation 


d  = 


p 


Ml 


d  T  Up 


where 


\\d\\l=  {xMyM 

and  vp  represents  the  sum  of  the  perturbation  accelerations.  Let  the  state  (s)  be 


Si  = 


then 


Si 


X* 

Zi 

Xi 

Yi 


Zi 


0 

0 

0 


0 

0 

0 

0 


0 

0 


Xi 


Yi 


Z^ 

Xi 

Yi 

Zi 


0 

0 

0 

0 

0 


IMilli 


1  0  0 
0  1  0 
0  0  1 
0  0  0 
0  0  0 
0  0  0 


X* 

0 

0 

Xi 

+ 

0 

^Px 

Yi 

^Py 

Zr_ 

^Pz 
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0.2  The  J2  Perturbation 

Because  gravity  is  a  conservative  force,  it  can  be  derived  from  the  gradient  of  a 
scalar  potential  function  [51].  Starting  with  Equation  (55),  Appendix  C  without  any 
control  inputs 

^  —pd 

d  =  — - h  Qjjj 

IMIli 
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Expanding  in  the  inertial  reference  frame  shown  in  Figure  84,  Appendix  C,  and  noting 


that  \d\  =  d=  \JX’^  +  where  X,Y ,  and  Z  are  the  position  of  the  satellite  in 

the  inertial  frame  I JK, 


d  = 


Noting  that 


-nX 


+ 


[X2  +  F2  +  Z2] 
—liZ 


j  +  d'pi 

2 


-IjlY 


[A2  +  Y^  +  Z2] 


T  + 
2 


J 


[X2  +  y2  +  ^2] 


J  +  dpK 
2 


K 


\/A:2  +  F2  +  ^2  ^  (0)  -  ^[l/2(\/A2  +  y2  +  Z2)-V2]2x' 


-^x 


[A2  +  y2  +  ^2] 

where  5  is  a  potential  function  such  that 


A2  +  y2  +  Z2 
J  +  ^pl 


Ojpj 


XB  =  dp 


Similarly 


d 

W 

A 


I  +  - 


-^Y 


[x^  +  Y^  +  z^y. 

—fiZ 

[X2  +  F2  +  ^2]| 


+ 


+  dpK 


Therefore  the  inertial  acceleration  is  the  gradient  of  the  potential  function  ^YB 


d  =  X[^YB 


The  B  term  comes  about  by  modeling  the  Earth  not  as  a  point  mass  but  as  an 
oblate  body  with  nonhomogeneous  mass  distribution  (Figure  95)  and  is  the  sum  of 
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Figure  95:  Oblate  Earth  Coordinate  Frame 


the  infinite  series 


i=2 


d 


JiPiisin  0)  +  XI  (  ^  if>)Pij{sm  0) 

i=i  ^  ^ 


(147) 


9?  =  jX  +  LVRe 


where  is  the  mean  eqnatorial  radius  of  the  Earth,  cj)  is  the  geocentric  latitnde  of  the 
satellite  (measured  from  the  equator) ,  A  is  the  geographical  longitude  (measured  from 
the  prime  meridian),  d  is  the  satellite  position  vector  magnitude,  cOe  is  the  rotation 
rate  of  the  Earth,  is  the  time  since  the  /  direction  lined  up  with  the  Greenwich 
meridian,  is  the  zonal  harmonic  coefficients  of  order  0,  is  a  Legendre  polynomial 
of  degree  n  and  order  0,  Pnm  is  a  Legendre  polynomial  of  degree  n  and  order  m,  Cnm 
is  the  tesseral  harmonic  coefficient  for  n  7^  m,  and  Snm  is  the  sectorical  harmonic 
coefficients  for  n  =  m. 

Measnrements  of  the  zonal,  tesseral  and  sectorial  coefficients  ( J„,  Cnm,  and  Snm) 
reveal  that  the  J2  term  is  at  least  400  times  larger  than  the  next  most  significant  term. 
Thus,  for  applications  such  as  satellite  reconfigurations  which  occur  over  relatively 
short  time  periods,  all  higher  terms  can  be  ignored.  Using  this  assumption,  Eqnation 
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(147)  reduces  to 


n=2 


R 

)  J„P„(sin0) 


=  ^  )  J2-P2(sm(/))  (148) 


where  J2  =  0.0010826  and  the  2nd  Legendre  polynomial  has  the  form  [24] 


ft(A')  =  -(3A^  -  1) 


Using  this  Legendre  polynomial  in  Equation  (148) 


Defining  the  constant 


^J2  —  7^J2R. 


2 

2-'--e 


simplifies  B  to 


^  =  ^J2  (  -^  [3(sin(/))2  -  1] 


By  geometry  in  Figure  95 


sin®  = 


d 


which  means 


'  z‘^ 

'if-' 


=  Mj2 


-3^2  +  (B 
'¥ 


Now  that  B  is  in  cartesian  coordinates,  dj^  is  found  by  taking  the  gradient  of  B 


dj„  —  VB  — 


—  R 

dx^ 

rfio 

A.  R 

—  /^^J2 

d^ (2d^)  —  (—3z‘^+d‘^)5d^  ^ 

1 - 

to  t 

rfio 

d5(-6Z+2df)-(-3Z2+d2)5^4| 

L  rfio 
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Simplifying 


^J2  —  /^^J2 


■  IbZ^X 

3A  ■ 

■  IbZ^x 

ZX  ' 

d\ 

“  If 

d\ 

~  If 

i5Z2y 

3y 

PJ2RI 

IhZ^Y 

AY 

d\ 

2 

df 

~  If 

15Z^ 

9Z 

15Z^ 

9Z 

d7 

~lf  . 

[  d7 

~lf  . 

(149) 


0.3  Atmospheric  Drag 

Atmospheric  drag  is  the  second  largest  perturbation  next  to  J2  for  satellites  in 
near-Earth  orbits  [59].  Unfortunately,  its  effects  are  considerably  harder  to  model 
accurately.  Unknowns  in  atmospheric  density,  spacecraft  ballistic  coefficient,  and 
relative  velocity  with  respect  to  the  atmosphere  all  contribute  to  errors  in  the  model. 
However,  a  simple  model  will  be  employed  in  order  to  give  a  sense  of  its  effects 
to  hovering  orbits.  The  standard  equation  for  acceleration  due  to  air  drag  is  well 
documented  in  the  literature  [26,50,59] 


*^drag 


ICpA 
2  m 


PdV^ 


Kel 


\\v. 


rel||2 


ICpA 
2  m 


where  Cp  is  the  coefficient  of  drag  (and  unitless),  A  is  the  area,  m  is  the  mass,  pd  is 
the  atmospheric  density,  and  Uei  is  the  satellite’s  velocity  relative  to  the  atmosphere 
and  equal  to 

brel  ^  ^atmos 


Due  to  the  Earth’s  rotation,  the  atmosphere  has  a  mean  velocity  of  [59] 


^atmos  ^  Af 


where  M  is  the  inertial  position  vector  of  the  satellite  and  the  rotation  of  the  Earth 


is  =  7.729  *  10 


-5 


rad 


.  Higher  order  effects  like  local  wind  and  density  variations 


will  be  ignored  in  this  simplified  model. 
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The  coefficient  of  drag,  area,  and  mass  are  usually  grouped  into  a  single  term 
called  the  ballistic  coefficient  {Be) 


Bc  = 


m 

CdA 


and  typically  have  a  range  of  [26] 


0.515 


m2 


<  Be  <  437 


m2 


The  last  piece  of  the  puzzle  is  a  model  of  atmospheric  density.  Although  a  number  of 


high  fidelity  stochastic  models  are  available,  we  can  get  the  accuracy  we  desire  with 


the  following  exponential  model  [59] 


hp  —  hd 

Pd  =  Poe  « 


where  po  is  the  nominal  density,  ho  is  the  reference  altitude,  hd  is  the  deputy’s  altitude, 
and  H  is  the  scale  height.  The  ho  and  H  values  are  determined  by  observations  and  are 
listed  in  Table  2  (verbatim  from  [59]).  This  model  uses  the  U.S.  Standard  Atmosphere 
(1976)  for  0  km,  CIRA-72  for  25-500  km,  and  CIRA-72  with  exospheric  temperature. 
Too  =  1000  K  for  500-1000  km.  The  atmospheric  drag  model  as  implemented  in  code 
is 


^  Rrel 

®drag  ot  D  P^e  ^  Irel 


2B, 


c 


(150) 
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Table  2:  Atmospheric  Model 


Altitude 

hd 

(km) 

Reference 
Altitude 
ho  (km) 

Nominal 
Density 
po  (kg/m^) 

Scale 
Height 
H  (km) 

Altitude 

hd 

(km) 

Reference 
Altitude 
ho  (km) 

Nominal 

Density 

Po  (kg/m2) 

Scale 
Height 
H  (km) 

0- 

25 

0 

1.225 

7.249 

150- 

180 

150 

2.070  X  10-9 

22.523 

25- 

30 

25 

3.899  X  10-2 

6.349 

180- 

200 

180 

5.464  X  10-^9 

29.740 

30- 

40 

30 

1.774  X  10-2 

6.682 

200- 

250 

200 

2.789  X  10-^9 

37.105 

40- 

50 

40 

3.972  X  10-^ 

7.554 

250- 

300 

250 

7.248  X  10-11 

45.546 

50- 

60 

50 

1.057  X  10-^ 

8.382 

300- 

350 

300 

2.418  X  10-11 

53.628 

60- 

70 

60 

3.206  X  10-^ 

7.714 

350  - 

400 

350 

9.518  X  10-12 

53.298 

70- 

80 

70 

8.770  X  10-^ 

6.549 

400  - 

450 

400 

3.725  X  10-12 

58.515 

80- 

90 

80 

1.905  X  10-^ 

5.799 

450  - 

500 

450 

1.585  X  10-12 

60.828 

90- 

100 

90 

3.396  X  10-® 

5.382 

500- 

600 

500 

6.967  X  10-13 

63.822 

100- 

no 

100 

5.297  X  10-^ 

5.877 

600- 

700 

600 

1.454  X  10-13 

71.835 

110- 

120 

no 

9.661  X  10-® 

7.263 

700  - 

800 

700 

3.614  X  10-1^ 

88.667 

120- 

130 

120 

2.438  X  10-® 

9.473 

800  - 

900 

800 

1.170  X  10-1^ 

124.64 

130- 

140 

130 

8.484  X  10-9 

12.636 

900  - 

1000 

900 

5.245  X  10-13 

181.05 

140  - 

150 

140 

3.845  X  10-9 

16.149 

>  1000 

1000 

3.019  X  10-13 

268.00 
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Appendix  P.  The  Gradient  of  the  Cost  Function 

Starting  with  the  cost  function  developed  in  Section  3.3,  Equation  (20) 


J 


AEJ  L 


+  lA/f  +  AVi  +  ...  +  AVf  +  AVf 


Tio  +  ^2  3  +  T3  4  + ...  + 


Tp 


we  can  calculate  the  gradient  with  respect  to  the  optimization  variables.  Let  the  total 
hover  time  be  Tt 

k 


Tp  +  7fc,fc+i  —  Tt 

1=1 

and  recall  that  AE^  is  a  function  of  the  previous,  current,  and  next  thrust  locations 
and  the  previous  and  next  time-of-flight 


AE/ ^i+i,  Ti^i+i) 


then  the  partials  of  J  with  respect  to  is 


dJ 

dipi 


Tt 


AVS 


d{AV,^)  a(AEh 
d'lpi  '  d'lpi 


1 


T2 

It 


TtAVS 


a(AE2)  ^  ^(AE^') 


d'lpi 


d'lpi 


(151) 


dJ 

dxp2 


T 


^AVS 


djAVp)  ,  djAVi)  8(AEb 
d'ip2  '  d'ip2  '  d'ip2 


-0 


n 


TtAVS 


diAV,^)  ^  djAVi)  ^  djAVi) 


d'ip2 


d'ip2 


d'ip2 
(152) 


dJ 

d'lpi 


Tt 


AV^ 


dipi  '  d'lpi  '  dipi 


-0 


Ti 


TrAVi 


djAVi_,)  ^  djAVj)  ^  djAVi,,) 


d'lpi 


dpJi  d'lpi 

(153) 


dJ 

d'lpk 


Tt 


1 

At/2 


9Ayp-i)  ,  dAVp) 

dipk  dipk 


n 


-0 


djAVl,)  djAVj)  djAVj) 

d'lpk  d'lpk  d'lpk 
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dJ 

d'lpk+i 


Tt-^ 


a(Ay2) 

sV’fc+i 


-0 


1 


n 


TtAV^ 


djAV,^)  ^  djAV^) 

di)k+i  dxpk+1 


Note  that  how  AV-^  and  its  partial  are  calculated  depends  on  the  entry  condition. 
Let 

k 

Y,AVI  =  AVt. 


i=l 


the  partials  with  respect  to  the  times-of-flight  are 


Tt 


Ays 


dJ 


noting  that  JTt  = 


,  a(Ay|) 

9Ti,2  9Ti,2 


It 


TtAVS 


d{AV,^)  ^  djAVi)  _  j 


on, 2 


dT,^2 


(154)' 


AVt2 


dJ 


Ays 


'a(AyS)  8(Av4y' 


OT,, 


2+1 


n 


TtAV^ 


^  a(AK+) 


,2+1 


^^2,2+1 

(155) 


aj 


Tt 


AV^ 


a(Ay?) 


AV^ 


A1/t2(1) 


5T, 


k,k-\-l 


T| 


TtAV^ 


diAV,^ 


d{AV^) 


dTi 


k,k+l 


dTi 


-  J 


k,k+l 


The  Tf  term  is  nonzero  in  only  one  exit  condition  where  it  is  set  equal  to  Ti,2;  therefore 
we  do  not  need  a  partial  with  respect  to  Tp. 


P.  1  The  Gradient  of  AV'^ 

A  few  preliminaries  are  required  before  taking  the  derivative  of  AV^.  Examining 
the  equation  for  R,  (Equation  (126),  Appendix  K)  we  see  only  a  handful  of  repeated 
forms,  thus  finding  the  gradient  is  not  as  complex  as  it  first  seems.  Let 

^  =  8  —  3nT  sin(nT)  —  8cos{nT) 
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and  its  partial  derivative  with  respect  to  time-of-flight  is  then 


dT 


—3n^T  cos(nT)  —  3nsin(nT)  +  8nsin(nT)  =  —3n^T  cos(nT)  +  5nsin(nT) 


Then  the  following  partials  are  (C  =  cos(nT)  and  S  =  sin(nT)) 


4S-3nT 


dT 


(8  -  3nTS  -  8C){4nC  -  3n)  -  (4,S  -  3nT){-3n^TC  +  3nS) 

(8  -  3nTS  -  8Cy 

-44n  +  56nC  +  2An‘^TS  -  9n^T^C  -  12nC^ 

(8  -  3nTS  -  8(7)2 


dT 


(8  -  3nTS  -  8C){-2nS)  -  (-2  +  2C){-3n^TC  +  bnS) 
(8  -  3nTS  -  8(7)2 
Qn^T  —  QnS  +  6nSC  —  Qn^TC 
(8  -  3nTS  -  8Cy 


^  (f )  _  (8  -  3nTS  -  8C)nC  -  S(-3n^TC  +  5nS) 

dT  ~  (8  -  3nTS  -  8Cf 

—bn  +  8nC  —  3nC^ 

(8  -  3nTS  -  8Cy 


4S-3nTC 

i 


dT 


(8  -  3nTS  -  8C){nC  +  3n‘^TS)  -  (4,S  -  3nTC){-3n^TC  +  bnS) 

(8  -  3nTS  -  8Cy 
-20n  +  8nC  +  2An^TS  +  l2nC‘^  -  9n^T^ 

(8  -  3nTS  -  8(7)2 


Q  ^  -14+6nT5+14C 

&T 


(8  -  3nTS  -  8C){6n^TC  -  8nS)  -  (-14  +  6nTS  +  14C')(-3n2rC'  +  bnS) 

(8  -  3nTS  -  8Cy 
6nS  -  QnSC  +  Qn^TC  -  Qn^T 
(8  -  3nTS  -  8(7)2 
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altogether 


dR 


dR 

W+ 


-44n+56nC-  +24n'^T- 5“  -Qn^T^-  C“  -12nC^- 

(8-3nT-5--8C-)2 

—6n^T~  +(>nS~  —6nS~  C~  +(>n^T~  C~ 

(8-3nT-5'--8C-)2 

20n-8nC~  -24n2r-  5"  -12nC2-  +9n^T^- 

(8-3nT-5--8C-)2 

6n2T-  -6n5-  +6nS-  Q- -Qn^T-  Q- 

(8-3nT-5--8C-)2 

0 

0 


0 

0 

20n-8nC+  -24n2r+  S+  -12nC^+  +9n^T2+ 

(8-3nr+S'+-8C+)2 

-6n2r+  +6n5+  -6n5+  C+  +6n2r+C+ 

(8-3nT+5+-8C+)2 

-44n+56nC++24n2r+  S+-9n^T^+  C+  -12nC^+ 

(8-3nT+5+-8C+)2 

6n^T~^  —6nS~^  +6nS'^C~^  —6n^T~^C'^ 

{8-3nT+ S+ -8C+)2 


6ii?T  —6nS  +QnS  C  —6n?T  C 

(8-3nT-5--8C-)2 

— 5n+8nC~  — 3nC^~ 

(8-3nT-5'--8C-)2 

-6n5-  +6n5-  C" -&'n?T-  C"  +6n2T- 

(8-3nT-5--8C-)2 

5n— 8nC~+3nC2~ 

(8-3nT-5--8C-)2 

0 


0 

0 

0 


6n5+-6n5+C++6n2T+C+-6n^T+ 

(8-3nT+S'+-8C+)2 

5n— 87x0+  +3nC2+ 

(8-3nT+5+-8C+)2 

— 6n2r++6n5'+— 6nS'+C"*"+6n2T"*'C"*' 
(8-3nT+5+-8C+)2 

— Sn+SnC"*"  — 

(8-3nT+5'+-8C+)2 


56) 


57) 


if  we  want  to  express  these  partials  in  terms  of  T,  we  simply  multiply  both  sides  by 


dT  j{^)  _ 

2it 

1 

'bn 

1 

1 

n 

and  substitute  nT  =  2ttT  from  Equation  (2),  yielding 

'  -887r+1127rC-  +9&iT^f-  §- -72TT^f^-C- -24^6^- 

24TT^f--12TTS-+12TTS-C--24TT^f-C- 

(8-67rf-S--8C-)2 

(8-67rf-S--8C'-)2 

-247r2f  - +127rS-  -127rS-  C"  +247r2T-C'“ 

—  IOtt+IGttC” 

(8-67rT-5--8C-)2 

(8-67rf-5--8C-)2 

OR 

407r-167rC- -96n^f-  5“  -247rC2-  +727r3T2- 

- 12715-  +12715-  C- -247i2t-  C-  +247i2f  - 

(8-67rT-5--8C'-)2 

(8-67if-5--8C-)2 

df-  ~ 

247r2T--127r5-+127r5-C'--247r2f-C- 

1071—1671 O'- +671(5^- 

(8-67rf-S--8C'-)2 

(8-67if-5--8C'-)2 

0 

0 

0 

0 

(158) 
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0 


0 


dk 

df+ 


0 

407r-167rC+  -967r^f'+  S+  -247rC^+ +727r^f^+ 

(8-67rf+S+-8C+)2 

-247r2f++127r:g+-127r:g+C'++247r2f+C+ 

(8-67rf+S+-8C+)2 

-887r+1127rC'++967r2f'+S+-727r3f'2+C+-247rC2+ 

(8-67rT+5+-8C+)2 

247r2f'+  -127r5++127r5+C'+  -247r2f +(7+ 

(8-67rT+5+-8C+)2 


0 

127r,g'+-127r,g'+C++247r^'f’+C+-247r^'f’+ 

(8-67rf'+:§+-8(5+)2 
IOtt— IGttC"*"  +67rC'2+ 

(8-67rf+S+-8(7+)2 

-  247r2  7^+ + 1 27r:g+  -  127r:g+  (7+  +247r2 1^+  (7+ 

(8-67rf+5+-8C+)2 
—  107r+167rC+— 67rC2+ 

(8-67rf+5+-8C+)2 

(159) 


We  are  now  prepared  to  take  the  derivative  of  with  respect  to  the  optimization 


variables.  Applying  the  theorem  found  in  Appendix  A. 4  and,  taking  the  partial 


derivative  with  respect  to  the  previous  point, 


2n^ 


^i—l  Ui—l  Vi  ^i+1  Vi+l 


RB! 


dxj-i 

dyi-i 

0 

0 

0 

0 


(160) 


where  the  partials  of  the  positions  are  functions  of  the  lobe  shape.  Equations  (42) 
and  (43)  from  Appendix  B 


Xi 


Vi 


7  cos  a  sin  P  + 
7  sin  a  sin  P  + 


_ TxTy  cos  Pi _ 

^ r2  cos2  {pi  -ri)+T^  sin^  {pi  -  rj) 

_ TxTy  sin  Pi _ 

t2  cos2  {Pi  -ri)  +  Tx  sin^  {pi  -  rp 
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The  derivatives  with  respect  to  ip  are  Equations  (47)  and  (48)  from  Appendix  B 


dxj 

dipi 

dyi 

dipi 


Ti 


-Ti  sin  ipi  -  rf  ^  ^  ^  sm{2ipi  -  2r])  cos  ipi 

X  y 

'^x  — 

Ti  cos  Ipi  -  rf  sm{2ipi  -  2r})  sin  ipi 

^x^y 

cos2  {ipi  -r})  +T^  sin^  {ipi  -  y) 


Similarly  the  partial  derivative  with  respect  to  the  current  point  is 


0 

0 


d  mn 

dipi 


2n^ 


^i—1  Vi—l  Vi  ^z+1  Vi+l 


RB! 


dxj 

dipi 

dyi 

dipi 


0 

0 


(161) 


and  with  respect  to  the  next  point: 


0 


0 


d{^Vp) 

dipi+i 


2n^ 


^i—1  Vi—l  Vi  ^z+1  J/z+1 


RR' 


0 

0 


dipi+i 

dyi+i 

.dipi+i_ 


(162) 
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Applying  the  theorem  found  in  Appendix  A. 5,  the  partial  derivative  with  respect  to 
the  previous  time-of-flight  is 


=  2n^ 


^i—l  Vi—l  Ui  ^i+1  Vi+l 


R 


dR' 


Xi-l 

Vi-i 

Xi 

Vi 

^i+l 

Vi+l 


(163) 


and  the  partial  with  respect  to  the  next  time-of-flight 


a  (Ay, 2) 


OR 


=  2n^ 


^i—1  Vi—l  Vi  ^i+1  Vi+l 


R 


dR' 

W+ 


Xi-l 

Vi-i 

Xi 

Vi 

Xi+l 

Vi+l 


(164) 


Recall  that  the  partial  derivatives  of  the  cost  function  J  requires  the  summation 
of  two  or  three  partials  of  AR  with  respect  to  '0  or  T.  When  implementing  these  in 
an  algorithm  it  becomes  useful  to  construct  the  following  matrices.  The  required 
summations  are  then  simply  the  row  sums. 
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AV^ 

AVi 

AVi 

Ay42 

■  ^  ^  AVil, 

AVI 

AV^ 

i’l 

9(AV2) 

a{Avi) 

d'tpi 

0 

0 

0 

0 

0 

9(Ab?) 

dipi 

a(Av^) 

dipi 

^2 

9(AV2) 

d'ip2 

a{Avi) 

d'ip2 

a{Avi) 

d'ip2 

0 

0 

0 

0 

0 

0 

^3 

0 

d{AVi) 

d^3 

d(AVi) 

d^3 

diAVi) 

d^3 

0 

0 

0 

0 

0 

^4 

0 

0 

d(AVi) 

dip4 

d{AVi) 

dlpA 

d{AVi) 

dipA 

0 

0 

0 

0 

0 

0 

0 

d{AVi) 

dtp5 

d{AVi) 

d^5 

0 

0 

0 

0 

i’k-2 

0 

0 

0 

0 

0 

diAVA2) 

dipk-2 

d(AVA,) 

d'4>k-2 

0 

0 

i’k-l 

0 

0 

0 

0 

0 

d{AVi_2) 

dipk-i 

mvAx) 

dipk-i 

9(Ab?) 

dipk-i 

0 

'Ipk 

0 

0 

0 

0 

0 

0 

aipk 

d{AVi) 

dipk 

9(AV2) 

dtpk 

'Ipk+l 

0 

0 

0 

0 

0 

0 

0 

d(AV^) 

d'ipk+1 

d(AV^) 

d'ipk+1 

The  last  two  entries  in  the  row  are  nonzero  only  in  the  repeating  hovering  orbit 
exit  condition.  The  partials  with  respect  to  time-of-flight 


A1^22 

AVi 

Alf 

AVi 

■  ■  ■  Al4l, 

AVti 

AVi 

Al/| 

Ti,2 

a(Av^) 

aTi,2 

a{Avi) 

9Ti,2 

0 

0 

0 

0 

0 

0 

9(AV2) 

9Ti,2 

2^2,3 

0 

d(AVi) 

dT2,3 

d(AVi) 

aT2,3 

0 

0 

0 

0 

0 

0 

^3,4 

0 

0 

d(AV^) 

dT3,4 

diAVi) 

9T3.4 

0 

0 

0 

0 

0 

^4,5 

0 

0 

0 

d{AVi) 

9T4,5 

d(AVi) 

9T4,5 

0 

0 

0 

0 

^5,6 

0 

0 

0 

0 

d(AVi) 

an, 6 

0 

0 

0 

0 

Tk-2,k-l 

0 

0 

0 

0 

0 

9{AV^_^) 

dTk-2,k-l 

d(AV^_,) 

dTk-2,k-i 

0 

0 

Tk-l,k 

0 

0 

0 

0 

0 

0 

d{AV^_,) 

aTk-\,k 

9(Ay2) 

aTk-\,k 

0 

Tk,k+1 

0 

0 

0 

0 

0 

0 

0 

d(AV^) 

dTk,k-\-i 

9(AV2) 
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Appendix  Q.  Matlab®  Algorithm 

The  following  appendix  contains  pseudocode  for  the  Matlab®  routine  used  to 
produce  the  results  presented. 


1.  Description:  Finds  and  plots  the  solution  for  staying  the  defined  lobe  while 
minimizing  the  total  AV  used  per  unit  of  time. 

2.  Inputs 


(a)  Entry/Exit  Condition 

(b)  Lobe  Parameters:  a,  (3,  7,  h,  r^,  Ty,  g 

(c)  Chief  Orbit  Parameters  (for  elliptical  orbits  only) 

(d)  ijj  Initial  Conditions  (determines  number  of  legs) 

3.  Algorithm 


(a) 

(b) 

(c) 


(d) 

(e) 


(f) 

(g) 

(h) 

(i) 

(j) 


Load  constraint  data  if  available 

Calculate  mean  motion  and  period  of  chief  (for  elliptical  chiefs) 
Calculate  Z  information 

i-  ^min  tII  cos /5 1 1  H  ^max  tII  C0S/3||  “h  H 

ii.  Max  Z  period  =  =  dcos"^(5mm/5max) 

Calculate  lobe  parameters:  Xmin,  Vxmin 

Construct  the  initial  condition  for  the  nonlinear  programming  algo¬ 
rithm  (set  T  =  0.01  chief  orbit  fractions  for  all) 

Ti^2,  ^2,3,  •  •  •  Tk^k+l] 


Tt 

A 


Run  EMINCON 

Calculate  number  of  legs  {kA  in  the  Z  direction: 

Generate  thrust  locations  and  the  time-of-flight  between  them 
Propagate  the  trajectory  between  thrust  locations 
Plot  data 
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